StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructSupport.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructSupport.cxx,v 1.30 2012/11/16 21:27:23 prindle Exp $
4  *
5  * Author: Jeff Porter
6  *
7  **********************************************************************
8  *
9  * Description: Simple helper class for calculating
10  * delta-rho, delta-rho/rho, and delta-rho/sqrt(rho)
11  * plus some other goodies
12  *
13  ***********************************************************************/
14 #include "StEStructSupport.h"
15 #include "Stiostream.h"
16 #include <sstream>
17 #include "Stsstream.h"
18 #include "TMath.h"
19 #include "TFile.h"
20 #include "TH1.h"
21 #include "TH2D.h"
22 
23 
24 const char* _pair_typename[] = {"Sib","Mix"};
25 const char* _pair_chargename[] = {"pp","pm","mp","mm"};
26 const char* _pair_ptweight[] = {"Pr","Su"};
27 int _pair_typemax=2;
28 int _pair_chargemax=4;
29 int _pair_totalmax=_pair_typemax*_pair_chargemax;
30 
31 ClassImp(StEStructSupport)
32 
33 //---------------------------------------------------------
34 bool StEStructSupport::goodName(const char* name){
35 
36  TString testName(_pair_typename[0]);
37  testName+=_pair_chargename[0];
38  testName+=name;
39  TObject* obj = NULL;
40  mtf->GetObject(testName.Data(),obj);
41  if(!obj) return false;
42 
43  return true;
44 }
45 //---------------------------------------------------------
46 bool StEStructSupport::goodName_zBuf(const char* name, int zBin){
47 
48  TString testName(_pair_typename[0]);
49  testName+=_pair_chargename[0];
50  testName+=name;
51  testName+="_zBuf_"; testName += zBin;
52  TObject* obj=NULL;
53  mtf->GetObject(testName.Data(),obj);
54  if(!obj) return false;
55 
56  return true;
57 }
58 
59 //---------------------------------------------------------
60 char* StEStructSupport::getFrontName(int itype){
61 
62  /* for itype=0-7, returns Sibpp, Sibpm, Sibmp, Sibmm, Mixpp, Mixpm, Mixmp, Mixmm */
63 
64  if(mtmpString) delete [] mtmpString;
65  mtmpString=new char[256];
66  std::ostringstream ts;
67  int j,k;
68  if(itype<_pair_chargemax){
69  j=0; k=itype;
70  } else {
71  j=1; k=itype-_pair_chargemax;
72  }
73 
74  ts<<_pair_typename[j]<<_pair_chargename[k];
75  strcpy(mtmpString,(ts.str()).c_str());
76  return mtmpString;
77 }
78 
79 //---------------------------------------------------------
80 const char* StEStructSupport::getTypeName(int itype){
81  return _pair_typename[itype];
82 }
83 
84 //---------------------------------------------------------
85 const char* StEStructSupport::getChargeSignName(int ics){
86  return _pair_chargename[ics];
87 }
88 
89 //---------------------------------------------------------
90 //
91 // Now the real class work
92 //
93 //---------------------------------------------------------
94 
95 StEStructSupport::StEStructSupport(TFile* tf, int bgmode, float* npairs) : mtf(tf), mbgMode(bgmode), mtmpString(NULL){
96 
97  //
98  // npairs is a normalization for when one cuts on say two (many) different
99  // ytyt slices and wants to compare the amplitudes, the generic normalization
100  // of sum of rho = 1. isn't sufficient
101  //
102 
103  msilent=false;
104  mapplyDEtaFix=false; // must set explicitly now
105  mDoSymmetrize = false;
106  mPairNormalization = false;
107  mPairWeighting = true;
108  mIdenticalPair = true;
109  mYtYtNormalization = false;
110  mYtYtVolumeNormalization = false;
111 
112  // Scan file for number of z-buffer bins
113  getNZBins();
114 
115  if(npairs){
116  mnpairs = new float[8];
117  for(int i=0;i<8;i++)mnpairs[i]=npairs[i];
118  } else {
119  mnpairs = 0;
120  }
121 
122 }
123 
124 StEStructSupport::~StEStructSupport(){
125  if(mtmpString) delete [] mtmpString;
126  if(mnpairs) delete [] mnpairs;
127 };
128 
129 //---------------------------------------------------------
130 int StEStructSupport::getNZBins(){
131  TString hname("NEventsSib_zBuf_");
132  int iZBin;
133  for (iZBin=1;iZBin<99;iZBin++) {
134  TString zBinName(hname.Data());
135  zBinName += iZBin-1;
136  TH1* tmp = NULL;
137  mtf->GetObject(zBinName.Data(),tmp);
138  if (!tmp) {
139  break;
140  }
141  }
142  mNumZBins = iZBin-1;
143  return mNumZBins;
144 };
145 //---------------------------------------------------------
146 // We have been calculating \Delta\rho/rho in different multiplicity/z bins
147 // then averaging using the number of pairs as a weight.
148 // This heavily favors slightly larger multiplicities..
149 // These get*Number routines are so I can try multiplicity weighting.
150 //
151 // Note: We allocate space for four floats, fill in the numbers, then return a pointer to them.
152 // Valgrind claims they are being used un-initialized. I think that is an issue
153 // with Valgrind. Might be nice to change how we return the numbers to get
154 // rid of those warnings.
155 float *StEStructSupport::getCommonNumber(int zBin) {
156  float *number = new float[4];
157  TString hname("");
158  TH1 *hpt[4];
159 
160  hname = "meanPtPA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[0]);
161  hname = "meanPtMA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[1]);
162  hname = "meanPtPB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[2]);
163  hname = "meanPtMB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[3]);
164 
165  // Want pp, pm, mp, mm
166  number[0] = hpt[0]->Integral() + hpt[2]->Integral();
167  number[1] = hpt[0]->Integral() + hpt[3]->Integral();
168  number[2] = hpt[1]->Integral() + hpt[2]->Integral();
169  number[3] = hpt[1]->Integral() + hpt[3]->Integral();
170  return number;
171 };
172 //---------------------------------------------------------
173 float *StEStructSupport::getCommonPairs(int zBin) {
174  float *pairs = new float[4];
175  for (int i=0;i<4;i++) {
176  pairs[i] = 0;
177  TString hname("");
178  hname += "Sib";
179  hname += _pair_chargename[i];
180  hname += "NDEtaDPhi_zBuf_";
181  hname += zBin;
182  TH1 *tmp = NULL;
183  mtf->GetObject(hname.Data(),tmp);
184  if (tmp) {
185  pairs[i] = tmp->Integral();
186  } else {
187  pairs[i] = 0;
188  }
189  }
190  return pairs;
191 };
192 //---------------------------------------------------------
193 float *StEStructSupport::getChargeNumber(int zBin) {
194  float *number = new float[4];
195  TString hname("");
196  TH1 *hpt[4];
197 
198  hname = "meanPtPA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[0]);
199  hname = "meanPtMA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[1]);
200  hname = "meanPtPB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[2]);
201  hname = "meanPtMB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt[3]);
202 
203  // Want LS, US, CI, CD. Turns out to just be grand sum.
204  number[0] = hpt[0]->Integral() + hpt[2]->Integral();
205  number[0]+= hpt[0]->Integral() + hpt[3]->Integral();
206  number[0]+= hpt[1]->Integral() + hpt[2]->Integral();
207  number[0]+= hpt[1]->Integral() + hpt[3]->Integral();
208  number[1] = number[0];
209  number[2]+= number[0];
210  number[3] = number[0];
211  return number;
212 };
213 //---------------------------------------------------------
214 float *StEStructSupport::getChargePairs(int zBin) {
215  float *pairs = new float[4];
216  for (int i=0;i<4;i++) {
217  pairs[i] = 0;
218  TString hname("");
219  hname += "Sib";
220  hname += _pair_chargename[i];
221  hname += "NDEtaDPhi_zBuf_";
222  hname += zBin;
223  TH1 *tmp = NULL;
224  mtf->GetObject(hname.Data(),tmp);
225  if (tmp) {
226  pairs[i] = tmp->Integral();
227  } else {
228  pairs[i] = 0;
229  }
230  }
231  pairs[0] += pairs[3];
232  pairs[1] += pairs[2];
233  pairs[2] = pairs[0] + pairs[1];
234  pairs[3] = pairs[0] + pairs[1];
235  return pairs;
236 };
237 
238 //---------------------------------------------------------
239  /* 09/27/12 djp Now taking eta range from Cuts file. Had been interrogating
240  * binning class for detaMax. Extract from special purpose histogram instead.
241  */
242 double StEStructSupport::getCIdNdEtadPhi() {
243  // A == sum over zBins {nA+} + {nA-}
244  // B == sum over zBins {nB+} + {nB-}
245  // return 0.5 sqrt(A*B) / nEvents
246 
247  double retVal = 0;
248  TH1 *hNSum;
249  TString hSibName;
250 
251  double nEvents = 0;
252  double nTracksA = 0;
253  double nTracksB = 0;
254  TH1 *hEta;
255  TString hname;
256  for (int iz=0;iz<mNumZBins;iz++) {
257  hSibName = "NEventsSib_zBuf_"; hSibName += iz; mtf->GetObject(hSibName.Data(),hNSum); nEvents += hNSum->Integral();
258 
259  hname = "etaPA_zBuf_"; hname += iz; mtf->GetObject(hname.Data(),hEta); nTracksA += hEta->Integral();
260  hname = "etaMA_zBuf_"; hname += iz; mtf->GetObject(hname.Data(),hEta); nTracksA += hEta->Integral();
261  hname = "etaPB_zBuf_"; hname += iz; mtf->GetObject(hname.Data(),hEta); nTracksB += hEta->Integral();
262  hname = "etaMB_zBuf_"; hname += iz; mtf->GetObject(hname.Data(),hEta); nTracksB += hEta->Integral();
263  }
264 
265  // The 0.5 is because we integrate tracks over two units of rapidity.
266  TH2D *hEtaPhi;
267  mtf->GetObject("EtaPhiRange",hEtaPhi);
268  double etaMin = -1;
269  double etaMax = +1;
270  if (hEtaPhi) {
271  TAxis *x = hEtaPhi->GetXaxis();
272  etaMin = x->GetXmin();
273  etaMax = x->GetXmax();
274  }
275  retVal = sqrt( nTracksA * nTracksB ) / (2*3.1415926 * nEvents * (etaMax-etaMin) );
276  return retVal;
277 }
278 //---------------------------------------------------------
279 double *StEStructSupport::getd2NdEtadPhi(int zBin, bool include2s) {
280  double* retVal = new double[6];
281  // (nA+)*(nB+), (nA+)*(nB-), (nA-)*(nB+), (nA-)*(nB-), (nA+)+(nA-)+(nB+)+(nB-)
282  // where the () means {d^2() \over d\eta d\phi}
283  // Note that in cases where we don't distinguish a and b
284  // (nA+)*(nB-) and (nA-)*(nB+) are the same and retVal[2] is not necessary.
285 
286 
287  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin; TH1* hNSum; mtf->GetObject(hSibName.Data(),hNSum);
288  double nEvents = hNSum->Integral();
289  retVal[5] = nEvents;
290 
291  double nTracksAP, nTracksAM, nTracksBP, nTracksBM;
292  double dNdEtaAP, dNdEtaAM, dNdEtaBP, dNdEtaBM;
293  TH1 *hEta;
294  TString hname;
295 
296  hname = "etaPA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hEta); nTracksAP = hEta->Integral();
297  hname = "etaMA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hEta); nTracksAM = hEta->Integral();
298  hname = "etaPB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hEta); nTracksBP = hEta->Integral();
299  hname = "etaMB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hEta); nTracksBM = hEta->Integral();
300 
301  // The 0.5 is because we integrate tracks over two units of rapidity.
302  TH2D *hEtaPhi;
303  mtf->GetObject("EtaPhiRange",hEtaPhi);
304  double etaMin = -1;
305  double etaMax = +1;
306  if (hEtaPhi) {
307  TAxis *x = hEtaPhi->GetXaxis();
308  etaMin = x->GetXmin();
309  etaMax = x->GetXmax();
310  }
311 
312  dNdEtaAP = nTracksAP / (nEvents * (etaMax-etaMin) );
313  dNdEtaAM = nTracksAM / (nEvents * (etaMax-etaMin) );
314  dNdEtaBP = nTracksBP / (nEvents * (etaMax-etaMin) );
315  dNdEtaBM = nTracksBM / (nEvents * (etaMax-etaMin) );
316  retVal[0] = dNdEtaAP*dNdEtaBP / pow(2*3.1415926,2);
317  retVal[1] = dNdEtaAP*dNdEtaBM / pow(2*3.1415926,2);
318  retVal[2] = dNdEtaAM*dNdEtaBP / pow(2*3.1415926,2);
319  retVal[3] = dNdEtaAM*dNdEtaBM / pow(2*3.1415926,2);
320 
321  // There are factors of 2 when expanding the full correlation term that
322  // could be included in loops over pairs, but appear not to have been.
323  // Fix that up here.
324  if(include2s) {
325  if (mIdenticalPair) {
326  retVal[4] = (dNdEtaAP + dNdEtaAM) / (2*3.1415926);
327  } else {
328  retVal[0] *= 2;
329  retVal[2] *= 2;
330  retVal[3] *= 2;
331  retVal[4] = (dNdEtaAP + dNdEtaAM + dNdEtaBP + dNdEtaBM) / (2*3.1415926);
332  }
333  }
334 
335  return retVal;
336 }
337 //---------------------------------------------------------
338 // These are the values of d2N/dEtadPhi used when combining \Delta\rho/\rho_{ref}
339 // First four values are product A+B+, A+B-, A-B+, A-A-.
340 // Fifth value is A+ + A- + B+ + B- (or one have of that, look in getd2NdEtadPhi).
341 double *StEStructSupport::getScaleFactors() {
342  double* retVal = new double[6];
343  for (int iType=0;iType<6;iType++) {
344  retVal[iType] = 0;
345  }
346 
347  double *d2NdEtadPhi;
348  for (int iz=0;iz<mNumZBins;iz++) {
349  d2NdEtadPhi = getd2NdEtadPhi(iz);
350  for (int iType=0;iType<5;iType++) {
351  retVal[iType] += d2NdEtadPhi[iType] * d2NdEtadPhi[5];
352  }
353  retVal[5] += d2NdEtadPhi[5];
354  delete [] d2NdEtadPhi;
355  }
356  for (int iType=0;iType<5;iType++) {
357  retVal[iType] /= retVal[5];
358  }
359  return retVal;
360 }
361 double *StEStructSupport::getScaleFactors(int zBin) {
362  return getd2NdEtadPhi(zBin);
363 }
364 //---------------------------------------------------------
365 double *StEStructSupport::getptHat() {
366  double* retVal = new double[4];
367  double d2N_Tot = 0;
368  double *d2NdEtadPhi;
369  for (int iType=0;iType<4;iType++) {
370  retVal[iType] = 0;
371  }
372  double *ptHat;
373  for (int iz=0;iz<mNumZBins;iz++) {
374  ptHat = getptHat(iz);
375  d2NdEtadPhi = getd2NdEtadPhi(iz);
376  for (int iType=0;iType<4;iType++) {
377  retVal[iType] += ptHat[iType] * d2NdEtadPhi[5];
378  }
379  d2N_Tot += d2NdEtadPhi[5];
380  delete [] d2NdEtadPhi;
381  delete [] ptHat;
382  }
383  for (int iType=0;iType<4;iType++) {
384  retVal[iType] /= d2N_Tot;
385  }
386  return retVal;
387 }
388 //---------------------------------------------------------
389 double *StEStructSupport::getptHat(int zBin) {
390  double* retVal = new double[4];
391  TH1 *hpt;
392  TString hname;
393 
394  hname = "meanPtPA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt); retVal[0] = hpt->GetMean();
395  hname = "meanPtMA_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt); retVal[1] = hpt->GetMean();
396  hname = "meanPtPB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt); retVal[2] = hpt->GetMean();
397  hname = "meanPtMB_zBuf_"; hname += zBin; mtf->GetObject(hname.Data(),hpt); retVal[3] = hpt->GetMean();
398 
399  return retVal;
400 }
401 //---------------------------------------------------------
402 int StEStructSupport::histogramExists(const char* name, int zBin){
403  TH2D* tmpF;
404  TString hname(getFrontName(0)); hname += name; hname+= "_zBuf_"; hname+= zBin;
405  mtf->GetObject(hname.Data(),tmpF);
406  if (tmpF) {
407  return 1;
408  }
409  return 0;
410 }
411 TH2D** StEStructSupport::getHists(const char* name, int zBin){
412  TH2D** retVal=NULL;
413  // If zBin = 0 check for name without with zBuf first.
414  if ((0 == zBin) && !goodName(name)) {
415  if (!goodName_zBuf(name,zBin)) {
416  return retVal;
417  }
418  }
419 
420  retVal = new TH2D*[8];
421  for (int i=0;i<_pair_totalmax;i++) {
422  TString hname(getFrontName(i)); hname+=name;
423  mtf->GetObject(hname.Data(),retVal[i]);
424  if (!retVal[i]) {
425  TString hname(getFrontName(i)); hname += name; hname+= "_zBuf_"; hname+= zBin;
426  mtf->GetObject(hname.Data(),retVal[i]);
427  }
428  }
429  return retVal;
430 }
431 //---------------------------------------------------------
432 TH2D** StEStructSupport::getLocalClones(const char* name, int zBin){
433  // A note on Sumw2. It appears that by default root will use the square root of the
434  // bin contents as the error. Not correct when we scale bin contents.
435  // One can use Sumw2 to trigger error propogation. However, this seems
436  // to be done implicitly in StEStructHAdd so we don't need to do it again.
437 
438  // hset contains pointers to in-memory histograms, _not_ copies.
439  TH2D** hset=getHists(name,zBin);
440  if(!hset) return (TH2D**)NULL;
441 
442  // make local clones
443  TH2D** hlocal=new TH2D*[_pair_totalmax];
444  for(int i=0;i<_pair_totalmax;i++) {
445  hlocal[i]=(TH2D*)hset[i]->Clone();
446 // hlocal[i]->Sumw2(); // trigger error propogation
447  }
448  // Delete array containing pointers, but not objects being pointed to.
449  delete [] hset;
450 
451  return hlocal;
452 }
453 //-----------------------------------------------------
454 void StEStructSupport::rescale(TH2D** hists) {
455  // Divide hists by bin widths, divide by Nevents.
456  // Bin by bin errors should be scaled properly, but won't include errors on the scale factors.
457 
458  if(!mtf) return;
459 
460  //>>>>>Need to rethink how to handle old cases where _zBuf_n is not in histogram name!!!!!
461  TH1 *hNum;
462  double nSib = 0;
463  double nMix = 0;
464  for (int zBin=0;zBin<mNumZBins;zBin++) {
465  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin; mtf->GetObject(hSibName.Data(),hNum);
466  nSib += hNum->Integral();
467  TString hMixName("NEventsMix_zBuf_"); hMixName += zBin; mtf->GetObject(hMixName.Data(),hNum);
468  nMix += hNum->Integral();
469  }
470 
471  for (int i=0;i<4;i++) {
472 // cout << "\ni=" << i << " Integral: " << hists[i]->Integral();
473  if (hists[i]->Integral() > 0) {
474  // dividing by ex*ey converts all input hists from counts to densities, so all operations and final result should also be density.
475  double dx = (hists[i]->GetXaxis())->GetBinWidth(1);
476  double dy = (hists[i]->GetYaxis())->GetBinWidth(1);
477  double binFactor = dx*dy;
478  hists[i]->Scale(1.0/(nSib*binFactor));
479  if (i==0 && !msilent) {
480  cout << "Scaling with Nevents " << nSib << " and binFactor " << binFactor << endl;
481  }
482 
483  if (mPairNormalization) {
484  // This is original normalization. Average value of \rho_{sib}/\rho_{ref} should be one.
485  TH2D *tmp = (TH2D *) hists[i]->Clone();
486  tmp->Divide(hists[i+4]);
487  double scale = tmp->Integral() / (tmp->GetNbinsX() * tmp->GetNbinsY());
488  hists[i+4]->Scale(scale);
489  delete tmp;
490  } else if (mYtYtNormalization) {
491  // Scale so hists[i+4]->Integral() = hists[i]->Integral()
492  // Seems that bins with small number of counts (the case for large Yt in YtYt plots)
493  // can skew the scaling when we scale as above.
494  hists[i+4]->Scale(hists[i]->Integral()/hists[i+4]->Integral());
495  } else if (mYtYtVolumeNormalization) {
496  // Scale so integral of \Delta\rho/sqrt(\rho) = 0;
497  TH2D *tmpA = (TH2D *) hists[i]->Clone();
498  tmpA->Reset();
499  TH2D *tmpB = (TH2D *) hists[i]->Clone();
500  tmpB->Reset();
501  for (int ix=1;ix<=tmpA->GetNbinsX();ix++) {
502  for (int iy=1;iy<=tmpA->GetNbinsY();iy++) {
503  double valS = hists[i]->GetBinContent(ix,iy);
504  double valR = hists[i+4]->GetBinContent(ix,iy);
505  if (valR > 0) {
506  tmpA->SetBinContent(ix,iy,valS/sqrt(valR));
507  } else {
508  tmpA->SetBinContent(ix,iy,0);
509  }
510  tmpB->SetBinContent(ix,iy,sqrt(valR));
511  }
512  }
513  double num = tmpA->Integral();
514  double den = tmpB->Integral();
515  if (den > 0) {
516  hists[i+4]->Scale(num/den);
517  } else {
518  hists[i+4]->Scale(0);
519  }
520  delete tmpA;
521  delete tmpB;
522  } else {
523  // We know how many sibling and mixed events there were.
524  // In this normalization we should be able to integrate (weighted with an
525  // appropriate kernel) to get fluctuations at larger scales.
526  if (nMix > 0) {
527  // For sibling the number of pairs are;
528  // n(n-1)/2 for ++ and --
529  // n^+ * n^- for +- (and -+)
530  // For mixed the number of pairs are
531  // n^(1) * n^(2) for ++ and --
532  // n^(+1) * n^(-1) + n^(-1) * n^(+2) for +-
533  // Thus we have about twice as many pairs per mixed event as sibling event.
534  hists[i+4]->Scale(0.5/(nMix*binFactor));
535  } else {
536  hists[i+4]->Scale(0);
537  }
538  }
539  }
540  }
541 }
542 
543 //-----------------------------------------------------
544 void StEStructSupport::rescale(TH2D** hists, int zBin) {
545  // Divide hists by bin widths, divide by Nevents.
546  // Bin by bin errors should be scaled properly, but won't include errors on the scale factors.
547 
548  if(!mtf) return;
549 
550  //>>>>>Need to rethink how to handle old cases where _zBuf_n is not in histogram name!!!!!
551  TH1 *hNum;
552  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin; mtf->GetObject(hSibName.Data(),hNum);
553  double nSib = hNum->Integral();
554  TString hMixName("NEventsMix_zBuf_"); hMixName += zBin; mtf->GetObject(hMixName.Data(),hNum);
555  double nMix = hNum->Integral();
556 
557  for (int i=0;i<4;i++) {
558 // cout << "\ni=" << i << " Integral: " << hists[i]->Integral();
559  if (hists[i]->Integral() > 0) {
560  // dividing by ex*ey converts all input hists from counts to densities, so all operations and final result should also be density.
561  double dx = (hists[i]->GetXaxis())->GetBinWidth(1);
562  double dy = (hists[i]->GetYaxis())->GetBinWidth(1);
563  double binFactor = dx*dy;
564  hists[i]->Scale(1.0/(nSib*binFactor));
565  if (i==0 && !msilent) {
566  cout << "Scaling with Nevents " << nSib << " and binFactor " << binFactor << endl;
567  }
568 
569  if (mPairNormalization) {
570  // This is original normalization. Average value of \rho_{sib}/\rho_{ref} should be one.
571  TH2D *tmp = (TH2D *) hists[i]->Clone();
572  tmp->Divide(hists[i+4]);
573  double scale = tmp->Integral() / (tmp->GetNbinsX() * tmp->GetNbinsY());
574  hists[i+4]->Scale(scale);
575  delete tmp;
576  } else if (mYtYtNormalization) {
577  // Scale so hists[i+4]->Integral() = hists[i]->Integral()
578  // Seems that bins with small number of counts (the case for large Yt in YtYt plots)
579  // can skew the scaling when we scale as above.
580  hists[i+4]->Scale(hists[i]->Integral()/hists[i+4]->Integral());
581  } else if (mYtYtVolumeNormalization) {
582  // Scale so integral of \Delta\rho/sqrt(\rho) = 0;
583  TH2D *tmpA = (TH2D *) hists[i]->Clone();
584  tmpA->Reset();
585  TH2D *tmpB = (TH2D *) hists[i]->Clone();
586  tmpB->Reset();
587  for (int ix=1;ix<=tmpA->GetNbinsX();ix++) {
588  for (int iy=1;iy<=tmpA->GetNbinsY();iy++) {
589  double valS = hists[i]->GetBinContent(ix,iy);
590  double valR = hists[i+4]->GetBinContent(ix,iy);
591  if (valR > 0) {
592  tmpA->SetBinContent(ix,iy,valS/sqrt(valR));
593  } else {
594  tmpA->SetBinContent(ix,iy,0);
595  }
596  tmpB->SetBinContent(ix,iy,sqrt(valR));
597  }
598  }
599  double num = tmpA->Integral();
600  double den = tmpB->Integral();
601  if (den > 0) {
602  hists[i+4]->Scale(num/den);
603  } else {
604  hists[i+4]->Scale(0);
605  }
606  delete tmpA;
607  delete tmpB;
608  } else {
609  // We know how many sibling and mixed events there were.
610  // In this normalization we should be able to integrate (weighted with an
611  // appropriate kernel) to get fluctuations at larger scales.
612  if (nMix > 0) {
613  hists[i+4]->Scale(0.5/(nMix*binFactor));
614  } else {
615  hists[i+4]->Scale(0);
616  }
617  }
618  }
619  }
620 }
621 
622 //---------------------------------------------------------
623 TH2D** StEStructSupport::getPtHists(const char* name, int zBin){
624 
625  TH2D** retVal=NULL;
626  // If zBin = 0 check for name without with zBuf first.
627  // These is a possiblitiy that the mean pt histograms are not in the file
628  // although the number correlations are.
629  // Add "Pr" in front of name before checking for existance of histogram.
630  TString chkName("Pr"); chkName += name;
631  if((0 == zBin) && !goodName(chkName.Data())) {
632  if (!goodName_zBuf(chkName.Data(),zBin)) {
633  return retVal;
634  }
635  }
636 
637  retVal=new TH2D*[32];
638 
639  for(int i=0;i<_pair_totalmax;i++){
640  TString hname(getFrontName(i)),hprname(getFrontName(i)),hpaname(getFrontName(i)),hpbname(getFrontName(i));
641  hname+=name; hname+="_zBuf_"; hname+=zBin;
642  hprname+="Pr"; hprname+=name; hprname+="_zBuf_"; hprname+=zBin;
643  hpaname+="Pa"; hpaname+=name; hpaname+="_zBuf_"; hpaname+=zBin;
644  hpbname+="Pb"; hpbname+=name; hpbname+="_zBuf_"; hpbname+=zBin;
645  mtf->GetObject(hname.Data(),retVal[i]);
646  mtf->GetObject(hprname.Data(),retVal[i+8]);
647  mtf->GetObject(hpaname.Data(),retVal[i+16]);
648  mtf->GetObject(hpbname.Data(),retVal[i+24]);
649  // Adjust errors on Pr, Pa and Pb. Current error is sqrt(bin content).
650  // Error on p_t is roughly \sigma_{p_t} = 0.015 * p_t
651  // We kind of interchange sums of squares and squares of sums to get following.
652  double n, pa, pb, er, ea, eb;
653  for (int ix=1;ix<=retVal[i+8]->GetNbinsX();ix++) {
654  for (int iy=1;iy<=retVal[i+8]->GetNbinsY();iy++) {
655  // Adjust sibling histograms for p_t uncertainties
656  // First four groups of retVal are ++, +-, -+ and -- for sibling.
657  // Second four are for mixed.
658  n = retVal[i]->GetBinContent(ix,iy);
659  pa = retVal[i+16]->GetBinContent(ix,iy);
660  pb = retVal[i+24]->GetBinContent(ix,iy);
661  er = (0.015*pa*pb/n)*sqrt(2/n);
662  ea = (0.015*pa)/sqrt(n);
663  eb = (0.015*pb)/sqrt(n);
664  retVal[i+8]->SetBinError(ix,iy,er);
665  retVal[i+16]->SetBinError(ix,iy,ea);
666  retVal[i+24]->SetBinError(ix,iy,eb);
667  }
668  }
669  }
670 
671  return retVal;
672 }
673 //---------------------------------------------------------------
674 TH2D** StEStructSupport::getPtClones(const char* name, int zBin){
675  // A note on Sumw2. See note in getClones above. I don't think we have done the
676  // >>>>> errors properly in StEStructHAdd yet. Need to look into this.
677 
678  // hset contains pointers to in-memory histograms, _not_ copies.
679  TH2D** hset=getPtHists(name,zBin);
680  if(!hset) return (TH2D**)NULL;
681 
682  // make local clones
683  TH2D** hlocal=new TH2D*[_pair_totalmax*4];
684  for(int i=0;i<_pair_totalmax*4;i++) {
685  hlocal[i]=(TH2D*)hset[i]->Clone();
686 // hlocal[i]->Sumw2(); // trigger error propogation
687  }
688  // Delete array containing pointers, but not objects being pointed to.
689  delete [] hset;
690 
691  return hlocal;
692 }
693 //-----------------------------------------------------
694 void StEStructSupport::rescalePt(TH2D** hists, int zBin) {
695  // Divide hists by bin widths, divide by Nevents.
696 
697  if(!mtf) return;
698 
699  TH1 *hNum;
700  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin; mtf->GetObject(hSibName.Data(),hNum);
701  double nSib = hNum->GetEntries();
702  TString hMixName("NEventsMix_zBuf_"); hMixName += zBin; mtf->GetObject(hMixName.Data(),hNum);
703  double nMix = hNum->GetEntries();
704 
705  for(int i=0;i<4;i++) {
706  double dx = (hists[i]->GetXaxis())->GetBinWidth(1);
707  double dy = (hists[i]->GetYaxis())->GetBinWidth(1);
708  double binFactor = dx*dy;
709  double nSibPairs = hists[i]->Integral();
710  double nMixPairs = hists[i+4]->Integral();
711  double norm;
712  if (nMixPairs > 0) {
713  norm = nSibPairs/(nMixPairs*nSib*binFactor);
714  } else {
715  norm = 0;
716  }
717  for (int j=0;j<4;j++) {
718  if (hists[i+4*j]->Integral() > 0) {
719  // divinding by ex*ey converts all input hists from counts to densities, so all operations and final result should also be density.
720  if(i==0 && !msilent) cout << "Scaling with Nevents " << nSib << " and binFactor " << binFactor << endl;
721  hists[i+8*j]->Scale(1.0/(nSib*binFactor));
722 
723  if (mPairNormalization) {
724  if (nMixPairs > 0) {
725  hists[i+4+8*j]->Scale(norm);
726  } else {
727  hists[i+4+8*j]->Scale(0);
728  }
729  } else {
730  if (nMix > 0) {
731  hists[i+4+8*j]->Scale(0.5/(nMix*binFactor));
732  } else {
733  hists[i+4+8*j]->Scale(0);
734  }
735  }
736  }
737  }
738  }
739 };
740 
741 //-----------------------------------------------------
742 void StEStructSupport::setSymmetrizeUS(bool symm) {
743  mDoSymmetrize = symm;
744 }
745 //-----------------------------------------------------
746 void StEStructSupport::symmetrizeUS(const char *name, TH2D** histos) {
747  // For pid mode the US X vs X and LS but different species histograms are not symmetrized.
748  // To form CD and CI combos we need the US to be symmetrized in the cases the LS are.
749 
750  // Histograms to be symmetrized are:
751  const char *symHistos[] = {"YtYt", "NYtYt", "PtPt",
752  "PhiPhi", "NPhiPhi", "PrPhiPhi", "PaPhiPhi", "PbPhiPhi",
753  "EtaEta", "PrEtaEta", "PaEtaEta", "PbEtaEta",
754  "EtaEtaSS", "PrEtaEtaSS", "PaEtaEtaSS", "PbEtaEtaSS",
755  "EtaEtaAS", "PrEtaEtaAS", "PaEtaEtaAS", "PbEtaEtaAS"};
756  // eight input histograms ++,+-,-+,-- for Sib and Mix
757  int symInt[] = {1,2, 5, 6};
758  for (int xy=0;xy<20;xy++) {
759  if (!strcmp(name,symHistos[xy])) {
760  for (int ih=0;ih<4;ih++) {
761  for (int ix=1;ix<=histos[symInt[ih]]->GetNbinsX();ix++) {
762  for (int iy=ix;iy<=histos[symInt[ih]]->GetNbinsY();iy++) {
763  double sum = histos[symInt[ih]]->GetBinContent(ix,iy) + histos[symInt[ih]]->GetBinContent(iy,ix);
764  double err = sqrt(pow(histos[symInt[ih]]->GetBinError(ix,iy),2) + pow(histos[symInt[ih]]->GetBinError(iy,ix),2));
765  histos[symInt[ih]]->SetBinContent(ix,iy,sum);
766  histos[symInt[ih]]->SetBinContent(iy,ix,sum);
767  histos[symInt[ih]]->SetBinError(ix,iy,err);
768  histos[symInt[ih]]->SetBinError(iy,ix,err);
769  }
770  }
771  }
772  }
773  }
774 };
775 //-----------------------------------------------------
776 void StEStructSupport::symmetrizePtUS(const char *name, TH2D** histos) {
777  // For pid mode the US X vs X and LS but different species histograms are not symmetrized.
778  // To form CD and CI combos we need the US to be symmetrized in the cases the LS are.
779 
780  // Histograms to be symmetrized are:
781  const char *symHistos[] = {"YtYt", "NYtYt", "PtPt",
782  "PhiPhi", "NPhiPhi", "PrPhiPhi", "PaPhiPhi", "PbPhiPhi",
783  "EtaEta", "PrEtaEta", "PaEtaEta", "PbEtaEta",
784  "EtaEtaSS", "PrEtaEtaSS", "PaEtaEtaSS", "PbEtaEtaSS",
785  "EtaEtaAS", "PrEtaEtaAS", "PaEtaEtaAS", "PbEtaEtaAS"};
786  int symInt[] = {1,2, 5, 6};
787  // 4 groups of 8 (Sibpp,Sibpm,Sibmp,Sibmm,Mixpp,Mixpm,Mixmp,Mixmm)
788  // 1st 8 are number, 2nd 8 are pt1*pt2, 3rd 8 are pt1 and 4th 8 are pt2
789  for (int xy=0;xy<20;xy++) {
790  if (!strcmp(name,symHistos[xy])) {
791  for (int ig=0;ig<32;ig+=8) {
792  for (int ih=0;ih<4;ih++) {
793  int histNum = ig + symInt[ih];
794  for (int ix=1;ix<=histos[histNum]->GetNbinsX();ix++) {
795  for (int iy=ix;iy<=histos[histNum]->GetNbinsY();iy++) {
796  double sum = histos[histNum]->GetBinContent(ix,iy) + histos[histNum]->GetBinContent(iy,ix);
797  double err = sqrt(pow(histos[histNum]->GetBinContent(ix,iy),2) + pow(histos[histNum]->GetBinContent(iy,ix),2));
798  histos[histNum]->SetBinContent(ix,iy,sum);
799  histos[histNum]->SetBinContent(iy,ix,sum);
800  histos[histNum]->SetBinError(ix,iy,err);
801  histos[histNum]->SetBinError(iy,ix,err);
802  }
803  }
804  }
805  }
806  }
807  }
808 };
809 
810 //-----------------------------------------------------
811 float* StEStructSupport::getNorms(TH2D** histArray){
812 
813  /* not really used any more but keep */
814 
815  float* retVal=NULL;
816  if(!histArray) return retVal;
817 
818  retVal=new float[8];
819  for(int i=0;i<_pair_totalmax; i++) retVal[i]=histArray[i]->Integral();
820 
821  return retVal;
822 }
823 
824 //---------------------------------------------------------
825 TH2D** StEStructSupport::buildCommonRatios(const char* name, float* sf){
826  return buildCommon(name,0,sf);
827 }
828 //---------------------------------------------------------
829 TH2D** StEStructSupport::buildCommonCFunctions(const char* name, float* sf){
830  return buildCommon(name,1,sf);
831 }
832 //---------------------------------------------------------
833 TH2D** StEStructSupport::buildCommonRFunctions(const char* name, float* sf){
834  return buildCommon(name,2,sf);
835 }
836 
837 //---------------------------------------------------------
838 TH2D** StEStructSupport::buildCommonRatios(const char* name, float* sf, int zBin){
839  return buildCommon(name,0,sf,zBin);
840 }
841 //---------------------------------------------------------
842 TH2D** StEStructSupport::buildCommonCFunctions(const char* name, float* sf, int zBin){
843  return buildCommon(name,1,sf,zBin);
844 }
845 //---------------------------------------------------------
846 TH2D** StEStructSupport::buildCommonRFunctions(const char* name, float* sf, int zBin){
847  return buildCommon(name,2,sf,zBin);
848 }
849 //---------------------------------------------------------
850 // Note for all buildXXX functions:
851 // We have three basic quantities, \rho_{sib}, \rho_{ref} and {d^2N \over d\eta d\phi}
852 // We use {d^2N \over d\eta d\phi}^2 as a proxy for \rho_{ref} in cases where we want
853 // an acceptance/efficiency independent number. (Ok, there is still an average
854 // efficiency that we have to divide by, and we probably want it at \eta=0.)
855 // When combining z-bins and/or centralities we want to add \Delta\rho, but
856 // each bin has a different efficiency correction. To account for this we define
857 // \hat\Delta\rho = {d^2N \over d\eta d\phi}^2 * \Delta\rho / \rho_{ref}.
858 // This takes role of \Delta\rho when we combine multiple bins.
859 // Option 0) \Delta\rho/\rho_{ref} -> \Sum_1^n{\hat\Delta\rho} / {(1/n) \Sum_1^n{dN/d\eta}}^2
860 // 1) \Delta\rho -> \hat\Delta\rho
861 // 2) \Delta\rho/sqrt(\rho_{ref}) -> \Sum_1^n{\hat\Delta\rho} / {(1/n) \Sum_1^n{dN/d\eta}}
862 // For YtYt we cannot use d2NdEtadPhi as proxy for \rho_{ref}.
863 // opt= 3, 4 and 5 combine \Delta\rho and \rho_{ref} directly.
864 // 3 \Delta\rho / rho_ref
865 // 4 \Delta\rho
866 // 5 \Delta\rh0 / sqrt(rho_ref)
867 
868 // Common -> {"PP","PM","MP","MM"};
869 TH2D** StEStructSupport::buildCommon(const char* name, int opt, float* sf) {
870  TH2D** retVal= buildCommon(name, opt, sf, 0);
871  if (!retVal) {
872  return retVal;
873  }
874 
875  // Fix up names and titles of histograms we will return.
876  for (int iType=0;iType<4;iType++) {
877  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
878  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
879  }
880  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
881  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
882  }
883  }
884  const char* oldName[4] = {"Sibpp","Sibpm","Sibmp","Sibmm"};
885  const char* newName[4] = {"PP ","PM ","MP ","MM "};
886  const char* oldTitle[4] = {"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-"};
887  const char* newTitle[4] = {"PP : ++ ","PM : +- ","MP: -+ ","MM : -- "};
888  for(int i=0;i<4;i++) {
889  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
890  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
891  }
892  for(int i=0;i<_pair_totalmax;i++) {
893  retVal[i]->Reset();
894  }
895 
896  if (opt < 3) {
897  // In the final scaling to \Delta\rho/\sqrt(\rho_{ref}) or \Delta\rho/\rho_{ref}
898  // one might want to use a {d^2N\over d\eta d\phi} for that specific
899  // charge combination. I always use the value appropriate for CI.
900 
901  // Note on combining centralites:
902  // Here we are taking ratios of sib/mix to cancel acceptance/efficiency before averaging.
903  // Combining centralities in this way will not give the same result as using a larger
904  // centrality bin to begin with, but may remove some of the systematic problems.
905 
906  // We weight each zBin by the total number of sibling pairs (integrated over (\eta_\Delta,\phi_\Delta))
907  // Divided by the number of pairs summed over zBins.
908 
909  // Twist is that if a Sibling OR Mixed histogram has a zero the ratio of root histograms
910  // gives 0. Want to ignore these in the averaging.
911  // Use histograms for weights, setting individual bins to 0 or weight as described above.
912  TH2D*** wScale = new TH2D**[mNumZBins];
913  TH2D** wTotal = new TH2D*[4];
914  for (int iType=0;iType<4;iType++) {
915  TString wName("sumWeighting"); wName += iType;
916  wTotal[iType] = (TH2D *) retVal[iType]->Clone(wName.Data());
917  }
918  for (int iz=0;iz<mNumZBins;iz++) {
919  wScale[iz] = getLocalClones(name, iz);
920  double weight;
921  for (int iType=0;iType<4;iType++) {
922  weight = wScale[iz][iType]->Integral();
923  for (int ix=1;ix<=wScale[iz][iType]->GetNbinsX();ix++) {
924  for (int iy=1;iy<=wScale[iz][iType]->GetNbinsY();iy++) {
925  double wTot = wTotal[iType]->GetBinContent(ix,iy);
926  if ((wScale[iz][iType]->GetBinContent(ix,iy) == 0) ||
927  (wScale[iz][iType+4]->GetBinContent(ix,iy) == 0)) {
928  wScale[iz][iType]->SetBinContent(ix,iy,0);
929  } else {
930  wScale[iz][iType]->SetBinContent(ix,iy,weight);
931  wTotal[iType]->SetBinContent(ix,iy,wTot+weight);
932  }
933  }
934  }
935  }
936  }
937  for (int iz=0;iz<mNumZBins;iz++) {
938  for (int iType=0;iType<4;iType++) {
939  wScale[iz][iType]->Divide(wTotal[iType]);
940  }
941  }
942 
943  // Now sum \Delta\rho / \rho_{ref} weighting by number of pairs in each zBin.
944  for (int iz=0;iz<mNumZBins;iz++) {
945  TH2D** tmpVal= buildCommon(name, opt, sf, iz);
946  for (int iType=0;iType<4;iType++) {
947  tmpVal[iType]->Multiply(wScale[iz][iType]);
948  retVal[iType]->Add(tmpVal[iType]);
949  delete tmpVal[iType];
950  }
951  delete [] tmpVal;
952  }
953 
954  // Have \Delta\rho / \rho_{ref}.
955  // Scale by some power of d^2N/d\eta d\phi (always using CI for this to be able to compare different charge combinations correctly).
956  // I was normalizing these \Delta\rho / \rho_{ref} by 1/4 thinking that would correct for
957  // adding four normalized histograms together. Actually that was approximately the algebra for
958  // combining to get CI, really wanted things like (n+/Nch)^2 which is close to 1/4.
959  // Move that scaling to code where we calculate CI, LS, US and CD (unless someone is relying
960  // on the old normalization for PP, PM, MP and MM in which case we move it back here.)
961  double *scale = getScaleFactors();
962  if (!msilent) {
963  cout << "@@@@@ In buildCommon: using d^2N/detadphi = " << scale[4] << endl;
964  }
965  for (int iType=0;iType<4;iType++) {
966  if (0 == opt) {
967  // Do nothing
968  } else if (1 == opt) {
969  retVal[iType]->Scale(scale[iType]);
970  } else if (2 == opt) {
971  retVal[iType]->Scale(sqrt(scale[iType]));
972  }
973  }
974 
975  // Delete space needed for scaling factors.
976  // Seem to be only deleting first four histograms most times when
977  // we getLocalClones. That really returns 8 histograms.
978  for (int iType=0;iType<8;iType++) {
979  for (int iz=0;iz<mNumZBins;iz++) {
980  delete wScale[iz][iType];
981  }
982  }
983  for (int iType=0;iType<4;iType++) {
984  delete wTotal[iType];
985  }
986  delete [] wScale;
987  delete [] wTotal;
988  } else {
989  double wIdentical[4] = {1, 1, 1, 1};
990  if (mIdenticalPair) {
991  wIdentical[1] = 2;
992  wIdentical[2] = 0;
993  }
994  // YtYt correlations
995  // I don't see any evidence for YtYt depending on the z-vertex position of the event.
996  // (I think there is clear evidence that (\eta_\Delta,\phi_\Delta) depends on z-vertex)
997  // It seems that in the ratio \rho{sib} / \rho_{ref} we have problems with bins where
998  // \rho_{ref} is small.
999  // For now just add zBins together, then do calculation.
1000 
1001  // buildCommon would scale each zBin. May want to refactor its functionality.
1002  // For now copy the appropriate part here.
1003  for (int iz=0;iz<mNumZBins;iz++) {
1004  TH2D** tmpVal = getLocalClones(name, iz);
1005  if (!tmpVal) {
1006  continue;
1007  }
1008  for (int iType=0;iType<8;iType++) {
1009  retVal[iType]->Add(tmpVal[iType]);
1010  delete tmpVal[iType];
1011  }
1012  delete [] tmpVal;
1013  }
1014 
1015  if (mDoSymmetrize) {
1016  symmetrizeUS(name,retVal);
1017  }
1018 
1019  if (mnpairs) { // manual scaling
1020  for(int i=0;i<8;i++) {
1021  retVal[i]->Scale(1.0/mnpairs[i]);
1022  }
1023  } else {
1024  rescale(retVal);
1025  }
1026 
1027  if (strstr(name,"DEta") || strstr(name,"SEta")) {
1028  fixDEtaGeometric((TH2**)retVal,8); // does nothing unless mapplyDEtaFix is set
1029  }
1030 
1031  float scf1=0.;
1032  float scf2=0.;
1033  if (sf) {
1034  scf1=sf[0];
1035  scf2=sf[1];
1036  }
1037  // if requested, scale bg to require correlations>=0 where statistics are large
1038  // This is important for Yt-Yt correlations. When we return to study those we
1039  // better double check this stuff is reasonable.
1040  if (1==mbgMode) {
1041  scaleBackGround(retVal[0],retVal[4],scf1);
1042  scaleBackGround(retVal[1],retVal[5],scf2);
1043  scaleBackGround(retVal[2],retVal[6],scf2);
1044  scaleBackGround(retVal[3],retVal[7],scf1);
1045  }
1046 
1047 
1048  // Calculate quantity we were asked for.
1049  for (int iType=0;iType<4;iType++) {
1050  // If rho and rho_ref are empty we skip the calculation
1051  if ((0 == retVal[iType]->Integral()) && (0 == retVal[iType]->GetMaximum()) &&
1052  (0 == retVal[iType+4]->Integral()) && (0 == retVal[iType+4]->GetMaximum())) {
1053  continue;
1054  }
1055  retVal[iType]->Scale(wIdentical[iType]);
1056  retVal[iType+4]->Scale(wIdentical[iType]);
1057  if (3 == opt) { // \Delta\rho/\rho_ref, Calculate \rho_sib/\rho_ref. Subtract 1 later.
1058  bool message = false;
1059  retVal[iType]->Divide(retVal[iType+4]);
1060  for (int ix=1;ix<=retVal[iType]->GetNbinsX();ix++) {
1061  for (int iy=1;iy<=retVal[iType]->GetNbinsY();iy++) {
1062  double rhoSib = retVal[iType]->GetBinContent(ix,iy);
1063  if (0 == rhoSib) {
1064  message = true;
1065  } else {
1066  retVal[iType]->SetBinContent(ix,iy,rhoSib-1);
1067  }
1068  }
1069  }
1070  if (message) {
1071  cout << "!!!!! Had at least one empty bin in Pt histogram " << name;
1072  cout << " type " << iType << " option 3 !!!!!" << endl;
1073  }
1074  } else if (4 == opt) { // \rho_sib - \rho_ref
1075  bool message = false;
1076  for (int ix=1;ix<=retVal[iType]->GetNbinsX();ix++) {
1077  for (int iy=1;iy<=retVal[iType]->GetNbinsY();iy++) {
1078  if (0 == retVal[iType]->GetBinContent(ix,iy)) {
1079  message = true;
1080  }
1081  }
1082  }
1083  if (message) {
1084  cout << "!!!!! Had at least one empty bin in Pt histogram " << name;
1085  cout << " type " << iType << " option 4!!!!!" << endl;
1086  }
1087  retVal[iType]->Add(retVal[iType+4],-1);
1088  } else if (5 == opt) { // \Delta\rho/sqrt(\rho_ref). Need to calculate errors by hand.
1089  TH1 *hNum;
1090  double nSib = 0;
1091  for (int zBin=0;zBin<mNumZBins;zBin++) {
1092  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin; mtf->GetObject(hSibName.Data(),hNum);
1093  nSib += hNum->Integral();
1094  }
1095  for (int ix=1;ix<=retVal[iType]->GetNbinsX();ix++) {
1096  for (int iy=1;iy<=retVal[iType]->GetNbinsY();iy++) {
1097  double rhoSib = retVal[iType]->GetBinContent(ix,iy);
1098  double rhoRef = retVal[iType+4]->GetBinContent(ix,iy);
1099  double eRhoSib = retVal[iType]->GetBinError(ix,iy);
1100  double eRhoRef = retVal[iType+4]->GetBinError(ix,iy);
1101  if (rhoRef > 0) {
1102  double eVal2 = eRhoSib*eRhoSib/rhoRef +
1103  0.25*eRhoRef*eRhoRef*pow(rhoSib/rhoRef+1,2)/rhoRef;
1104  double val = (rhoSib-rhoRef)/sqrt(rhoRef);
1105  retVal[iType]->SetBinContent(ix,iy,val);
1106  retVal[iType]->SetBinError(ix,iy,sqrt(eVal2));
1107  } else {
1108  double eVal2 = eRhoSib*eRhoSib + 1;
1109  retVal[iType]->SetBinContent(ix,iy,0);
1110  // Not sure about this error.
1111  retVal[iType]->SetBinError(ix,iy,sqrt(eVal2/nSib));
1112  }
1113  }
1114  }
1115  }
1116  }
1117 
1118  }
1119  return retVal;
1120 }
1121 //---------------------------------------------------------
1122 // Note that I (djp) modified the analysis code so that we only fill the -+
1123 // histograms for identified, different particles. So we get -+ filled for
1124 // pi-K+ as an example but not for unidentified particles (as in mode 1
1125 // and mode 3) where we don't know what the particle types are.
1126 
1127 // Ignore opt here, only use it for over-loading distinction.
1128 // For opt < 3 return \Delta\rho/rho_{ref}
1129 // For larger opt return rho and rho_{ref}.
1130 TH2D** StEStructSupport::buildCommon(const char* name, int opt, float* sf, int zBin) {
1131 
1132  /* builds hist types = ++,+-,-+,-- */
1133 
1134 
1135  // eight input histograms ++,+-,-+,-- for Sib and Mix
1136  TH2D** retVal = getLocalClones(name, zBin);
1137  if (!retVal) {
1138  return retVal;
1139  }
1140 
1141  if (mDoSymmetrize) {
1142  symmetrizeUS(name,retVal);
1143  }
1144 
1145  if (mnpairs) { // manual scaling
1146  for(int i=0;i<8;i++) {
1147  retVal[i]->Scale(1.0/mnpairs[i]);
1148  }
1149  } else {
1150  rescale(retVal, zBin); // Divide out bin size, number of events to approximate density.
1151  // Optionally scale number of mix pairs to be same as sibling.
1152  }
1153  if (strstr(name,"DEta") || strstr(name,"SEta")) {
1154  fixDEtaGeometric((TH2**)retVal,8); // does nothing unless mapplyDEtaFix is set
1155  }
1156 
1157  float scf1=0.;
1158  float scf2=0.;
1159  if (sf) {
1160  scf1=sf[0];
1161  scf2=sf[1];
1162  }
1163  // if requested, scale bg to require correlations>=0 where statistics are large
1164  // This is important for Yt-Yt correlations. When we return to study those we
1165  // better double check this stuff is reasonable.
1166  if (1==mbgMode) {
1167  scaleBackGround(retVal[0],retVal[4],scf1);
1168  scaleBackGround(retVal[1],retVal[5],scf2);
1169  scaleBackGround(retVal[2],retVal[6],scf2);
1170  scaleBackGround(retVal[3],retVal[7],scf1);
1171  }
1172 
1173  if (opt < 3) {
1174  // \Delta\rho/\rho_{ref}
1175  // Actually calculate \rho/\rho_{ref} - 1 so errors are correct.
1176  // Don't know a convenient way to subtract one from every bin!!!!
1177  // Note that when a bin in retVal[i] or retVal[i+4] is 0 the ratio
1178  // will be set to 0 and then we subtract 1. Check for this and don't subtract the 1.
1179  // In parent method we hopefully weight this bin with 0 so don't include in an average.
1180  for (int i=0;i<4;i++) {
1181  // If rho and rho_ref are empty we skip the calculation
1182  if ((0 == retVal[i]->Integral()) && (0 == retVal[i]->GetMaximum()) &&
1183  (0 == retVal[i+4]->Integral()) && (0 == retVal[i+4]->GetMaximum())) {
1184  continue;
1185  }
1186  retVal[i]->Divide(retVal[i+4]); // rho/rho_{ref}
1187  for (int ix=1;ix<=retVal[i]->GetNbinsX();ix++) {
1188  for (int iy=1;iy<=retVal[i]->GetNbinsY();iy++) {
1189  double val = retVal[i]->GetBinContent(ix,iy);
1190  if (val == 0) {
1191  retVal[i]->SetBinContent(ix,iy,0);
1192  } else {
1193  retVal[i]->SetBinContent(ix,iy,val-1); // delta-rho/rho_{ref}
1194  }
1195  }
1196  }
1197  }
1198  }
1199 
1200  return retVal;
1201 }
1202 //---------------------------------------------------------
1203 // We should not need to subtract mixed pairs in calculation of pt correlations.
1204 // It turns out to make a surprisingly big difference.
1205 
1206 // Combine zBuffers as we do for number correlations.
1207 
1208 // opt: 0 = delta-rho/rho_mix; 1 = delta-rho; 2 = delta-rho/sqrt(rho_mix);
1209 TH2D** StEStructSupport::buildPtCommon(const char* name, int opt, int subtract) {
1210  TH2D** retVal= buildPtCommon(name, opt, subtract, 0);
1211  if (!retVal) {
1212  return retVal;
1213  }
1214  // Fix up name and title of histograms we return
1215  for (int iType=0;iType<4;iType++) {
1216  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
1217  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
1218  }
1219  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
1220  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
1221  }
1222  }
1223  const char* oldName[4] = {"Sibpp","Sibpm","Sibmp","Sibmm"};
1224  const char* newName[4] = {"PP ","PM ","MP ","MM "};
1225  const char* oldTitle[4] = {"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-"};
1226  const char* newTitle[4] = {"PP : ++ ","PM : +- ","MP: -+ ","MM : -- "};
1227  for (int i=0;i<4;i++) {
1228  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
1229  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1230  retVal[i]->Reset();
1231  }
1232 
1233  // In the final scaling to \Delta\rho/\sqrt(\rho_{ref}) or \Delta\rho/\rho_{ref}
1234  // one might want to use a {d^2N\over d\eta d\phi} for that specific
1235  // charge combination. I always use the value appropriate for CI.
1236 
1237  // Note on combining centralites:
1238  // Here we are taking ratios of sib/mix to cancel acceptance/efficiency before averaging.
1239  // Combining centralities in this way will not give the same result as using a larger
1240  // centrality bin to begin with, but may remove some of the systematic problems.
1241 
1242  // We weight each zBin by the total number of sibling pairs (integrated over (\eta_\Delta,\phi_\Delta)).
1243  // First sum over zBins to get total number of pairs.
1244  TH2D** temp;
1245  double *wScale[4];
1246  wScale[0] = new double[mNumZBins];
1247  wScale[1] = new double[mNumZBins];
1248  wScale[2] = new double[mNumZBins];
1249  wScale[3] = new double[mNumZBins];
1250  double wTotal[4] = {0, 0, 0, 0};
1251  for (int iz=0;iz<mNumZBins;iz++) {
1252  temp = getLocalClones(name, iz);
1253  for (int iType=0;iType<4;iType++) {
1254  wScale[iType][iz] = temp[iType]->Integral();
1255  wTotal[iType] += wScale[iType][iz];
1256  delete temp[iType];
1257  }
1258  for (int iType=4;iType<8;iType++) {
1259  delete temp[iType];
1260  }
1261  delete [] temp;
1262  }
1263 
1264  // Now sum \Delta\rho / \rho_{ref} weighting by number of pairs in each zBin.
1265  for (int iz=0;iz<mNumZBins;iz++) {
1266  TH2D** tmpVal= buildPtCommon(name, opt, subtract, iz);
1267  for (int iType=0;iType<4;iType++) {
1268  if (0 != wTotal[iType]) {
1269  retVal[iType]->Add(tmpVal[iType], wScale[iType][iz]/wTotal[iType]);
1270  }
1271  delete tmpVal[iType];
1272  }
1273  delete [] tmpVal;
1274  }
1275 
1276  // Have weighted \Delta\rho / \rho_{ref}.
1277  // Scale by some power of d^2N/d\eta d\phi (always using CI for this to be able to compare different charge combinations correctly).
1278  double *scale = getScaleFactors();
1279  if (!msilent) {
1280  cout << "@@@@@ In buildPtCommon: using d^2N/detadphi = " << scale[4] << endl;
1281  }
1282  for (int iType=0;iType<4;iType++) {
1283  if (0 == opt) {
1284  // Do nothing
1285  } else if (1 == opt) {
1286  retVal[iType]->Scale(scale[iType]);
1287  } else if (2 == opt) {
1288  retVal[iType]->Scale(sqrt(scale[iType]));
1289  }
1290  }
1291 
1292  return retVal;
1293 }
1294 //---------------------------------------------------------
1295 TH2D** StEStructSupport::buildPtCommon(const char* name, int opt, int subtract, int zBin) {
1296 
1297  /* builds hist types = ++,+-,-+,-- */
1298 
1299  if(!mtf) return (TH2D**)NULL;
1300 
1301  // -- here we get 32 hists:
1302  // 4 groups of 8 (Sibpp,Sibpm,Sibmp,Sibmm,Mixpp,Mixpm,Mixmp,Mixmm)
1303  // 1st 8 are number, 2nd 8 are pt1*pt2, 3rd are pt1 and 4th are pt2
1304 
1305  TH2D** hlocal=getPtClones(name,zBin);
1306  if (!hlocal) {
1307  return hlocal;
1308  }
1309 
1310  if (mDoSymmetrize) {
1311  symmetrizePtUS(name,hlocal);
1312  }
1313  rescalePt(hlocal, zBin); // Divide out bin size, number of events to approximate density.
1314  // Optionally scale number of mix pairs to be same as sibling.
1315 
1316  // four returned hists
1317  TH2D** retVal= new TH2D*[4]; // 0=++ 1=+- 2=-+ 3=--
1318  TH2D** mixVal= new TH2D*[4]; // 0=++ 1=+- 2=-+ 3=--
1319  for (int i=0;i<4;i++) {
1320  retVal[i]=(TH2D*)hlocal[i]->Clone();// just get correct dimensions & names
1321  retVal[i]->Scale(0.); // zero the hists
1322  mixVal[i]=(TH2D*)hlocal[0]->Clone();
1323  mixVal[i]->Scale(0.); // zero the hists
1324  }
1325 
1326  if(strstr(name,"DEta") || strstr(name,"SEta")) {
1327  fixDEtaGeometric((TH2**)hlocal,32);
1328  }
1329  double *ptHat = getptHat(zBin);
1330  double ptHatA[4], ptHatB[4];
1331  ptHatA[0] = ptHat[0];
1332  ptHatA[1] = ptHat[0];
1333  ptHatA[2] = ptHat[1];
1334  ptHatA[3] = ptHat[1];
1335  ptHatB[0] = ptHat[2];
1336  ptHatB[1] = ptHat[2];
1337  ptHatB[2] = ptHat[3];
1338  ptHatB[3] = ptHat[3];
1339  delete [] ptHat;
1340  for (int i=0;i<4;i++) { // Note that i=0->3 is ++, +-, -+, --
1341  retVal[i]->Add(hlocal[i+ 8]);
1342  retVal[i]->Add(hlocal[i+16],-ptHatA[i]);
1343  retVal[i]->Add(hlocal[i+24],-ptHatB[i]);
1344  retVal[i]->Add(hlocal[i ],ptHatA[i]*ptHatB[i]);
1345 
1346  mixVal[i]->Add(hlocal[i+12]);
1347  mixVal[i]->Add(hlocal[i+20],-ptHatA[i]);
1348  mixVal[i]->Add(hlocal[i+28],-ptHatB[i]);
1349  mixVal[i]->Add(hlocal[i+ 4],ptHatA[i]*ptHatB[i]);
1350  }
1351 
1352  for (int i=0;i<4;i++) {
1353  if (subtract) {
1354  // Zero errors of mixed events. Rationale is that those errors are
1355  // already included in the sibling.
1356  for (int ix=1;ix<=mixVal[i]->GetNbinsX();ix++) {
1357  for (int iy=1;iy<=mixVal[i]->GetNbinsY();iy++) {
1358  mixVal[i]->SetBinError(ix,iy,0);
1359  }
1360  }
1361  retVal[i]->Add(mixVal[i],-1); // Subtract mixed event artifacts
1362  }
1363  retVal[i]->Divide(hlocal[i+4]);
1364  }
1365 
1366  // Free local histograms.
1367  for (int i=1;i<4;i++) {
1368  delete mixVal[i];
1369  }
1370  for(int i=0;i<_pair_totalmax*4;i++) {
1371  delete hlocal[i];
1372  }
1373  delete [] mixVal;
1374  delete [] hlocal;
1375 
1376  return retVal;
1377 }
1378 
1379 //---------------------------------------------------------
1380 TH2D** StEStructSupport::buildChargeTypeRatios(const char* name, float* sf){
1381  return buildChargeTypes(name,0,sf);
1382 }
1383 //---------------------------------------------------------
1384 TH2D** StEStructSupport::buildChargeTypeCFunctions(const char* name, float* sf){
1385  return buildChargeTypes(name,1,sf);
1386 }
1387 //---------------------------------------------------------
1388 TH2D** StEStructSupport::buildChargeTypeRFunctions(const char* name, float* sf){
1389  return buildChargeTypes(name,2,sf);
1390 }
1391 //---------------------------------------------------------
1392 // See comments before buildCommon about d2NdEtadPhi being turned into proxy for \rho_{ref}
1393 // ChargeTypes -> {"LS","US","CD","CI"}
1394 // For opt=0 calculate \Delta\rho / rho_ref
1395 // 1 \Delta\rho
1396 // 2 \Delta\rh0 / sqrt(rho_ref)
1397 // For YtYt we cannot use d2NdEtadPhi as proxy for \rho_{ref}.
1398 // opt= 3, 4 and 5 combine \Delta\rho and \rho_{ref} directly.
1399 // 3 \Delta\rho / rho_ref
1400 // 4 \Delta\rho
1401 // 5 \Delta\rh0 / sqrt(rho_ref)
1402 TH2D** StEStructSupport::buildChargeTypes(const char* name, int opt, float *sf) {
1403  // Scale ++, +-, -+ and -- by their respective {d^2N \over d\eta d\phi}^2 to create
1404  // \Delta\rho. Scale by appropriate amount of {d^2N \over d\eta d\phi} for opt.
1405  // Finally add together to form LS, US, CD and CI.
1406 
1407  // 1/14/08 djp This routine did exactly the same thing as buildCommon up to the point
1408  // of forming LS, US, CD and CI. Use that routine instead so we have fewer
1409  // places for bugs.
1410  TH2D** retVal = buildCommon(name, opt, sf);
1411  if (!retVal) {
1412  return retVal;
1413  }
1414 
1415  // Fix up names and titles of histograms we will return.
1416  for (int iType=0;iType<4;iType++) {
1417  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
1418  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
1419  }
1420  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
1421  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
1422  }
1423  }
1424  const char* oldName[4] = {"PP ","PM ","MP ","MM "};
1425  const char* newName[4] = {"LS ","US ","CD ","CI "};
1426  const char* oldTitle[4] = {"PP : ++ ","PM : +- ","MP: -+ ","MM : -- "};
1427  const char* newTitle[4] = {"LS : ++ + -- ","US : +- + -+ ","CD: ++ + -- - +- - -+ ","CI : ++ + -- + +- + -+ "};
1428  for(int i=0;i<4;i++) {
1429  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
1430  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1431  }
1432 
1433  // Now form LS, US, CD and CI
1434  double *scale = getScaleFactors();
1435  for (int iType=0;iType<4;iType++) {
1436  // Need to undo scaling with charge specific normalization and redo with CI
1437  if (0 == opt) {
1438  retVal[iType]->Scale(scale[iType]/pow(scale[4],2));
1439  } else if (1 == opt) {
1440  retVal[iType]->Scale(1.0);
1441  } else if (2 == opt) {
1442  retVal[iType]->Scale(sqrt(scale[iType])/scale[4]);
1443  }
1444  }
1445  delete [] scale;
1446  retVal[0]->Add(retVal[3]);
1447  if (mIdenticalPair) {
1448  retVal[1]->Scale(2);
1449  } else {
1450  retVal[1]->Add(retVal[2]);
1451  }
1452  retVal[2]->Add(retVal[0],retVal[1],1,-1);
1453  retVal[3]->Add(retVal[0],retVal[1]);
1454 
1455  return retVal;
1456 }
1457 //---------------------------------------------------------
1458 TH2D** StEStructSupport::buildChargeTypeRatios(const char* name, float* sf, int zBin){
1459  return buildChargeTypes(name,0,sf,zBin);
1460 }
1461 //---------------------------------------------------------
1462 TH2D** StEStructSupport::buildChargeTypeCFunctions(const char* name, float* sf, int zBin){
1463  return buildChargeTypes(name,1,sf,zBin);
1464 }
1465 //---------------------------------------------------------
1466 TH2D** StEStructSupport::buildChargeTypeRFunctions(const char* name, float* sf, int zBin){
1467  return buildChargeTypes(name,2,sf,zBin);
1468 }
1469 //---------------------------------------------------------
1470 // This is not used except by buildChargeTypeRatios, buildChargeTypeCFunctions and
1471 // buildChargeTypeRFunctions. This might be useful as a check on zBuffer dependence.
1472 TH2D** StEStructSupport::buildChargeTypes(const char* name, int opt, float *sf, int zBin) {
1473  // Scale ++, +-, -+ and -- by their respective {d^2N \over d\eta d\phi}^2 to create
1474  // \Delta\rho. Scale by appropriate amount of {d^2N \over d\eta d\phi} for opt.
1475  // Finally add together to form LS, US, CD and CI.
1476 
1477  TH2D** retVal= buildCommon(name, opt, sf, zBin);
1478  if (!retVal) {
1479  return retVal;
1480  }
1481 
1482  // Fix up names and titles of histograms we will return.
1483  for (int iType=0;iType<4;iType++) {
1484  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
1485  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
1486  }
1487  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
1488  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
1489  }
1490  }
1491  const char* oldName[4] = {"Sibpp","Sibpm","Sibmp","Sibmm"};
1492  const char* newName[4] = {"LS ","US ","CD ","CI "};
1493  const char* oldTitle[4] = {"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-"};
1494  const char* newTitle[4] = {"LS : ++ + -- ","US : +- + -+ ","CD: ++ + -- - +- - -+ ","CI : ++ + -- + +- + -+ "};
1495  for(int i=0;i<4;i++) {
1496  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
1497  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1498  }
1499 
1500  if (opt < 3) {
1501  // Have \Delta\rho/rho_{ref}.
1502  // Scale for desired option.
1503  double *d2NdEtadPhi = getd2NdEtadPhi(zBin);
1504  double sqrtRho = d2NdEtadPhi[4];
1505  delete [] d2NdEtadPhi;
1506  for (int iType=0;iType<4;iType++) {
1507  if (0 < sqrtRho) {
1508  if (1 == opt) {
1509  retVal[iType]->Scale(pow(sqrtRho,2));
1510  } else if (2 == opt) {
1511  retVal[iType]->Scale(sqrtRho);
1512  }
1513  } else {
1514  retVal[iType]->Scale(0.0);
1515  }
1516  }
1517  } else {
1518  // Scale according to opt.
1519  // Here we cannot use d2NdEtadPhi as proxy for rho_{ref}
1520  for (int iType=0;iType<4;iType++) {
1521  retVal[iType]->Add(retVal[iType+4],-1);
1522  for (int ix=1;ix<=retVal[iType]->GetNbinsX();ix++) {
1523  for (int iy=1;iy<=retVal[iType]->GetNbinsY();iy++) {
1524  double rhoRef = retVal[iType+4]->GetBinContent(ix,iy);
1525  if (1 < rhoRef) {
1526  double val = retVal[iType]->GetBinContent(ix,iy);
1527 // double err = retVal[iType]->GetBinError(ix,iy);
1528  if (3 == opt) {
1529  retVal[iType]->SetBinContent(ix,iy,val/rhoRef);
1530 // retVal[iType]->SetBinError(ix,iy,err/rhoRef);
1531  } else if (5 == opt) {
1532  retVal[iType]->SetBinContent(ix,iy,val/sqrt(rhoRef));
1533 // retVal[iType]->SetBinError(ix,iy,err/sqrt(rhoRef));
1534  }
1535  }
1536  }
1537  }
1538  }
1539  }
1540 
1541  // Now form LS, US, CD and CI
1542  retVal[0]->Add(retVal[3]);
1543  retVal[1]->Add(retVal[2]);
1544  retVal[2]->Add(retVal[0],retVal[1],1,-1);
1545  retVal[3]->Add(retVal[0],retVal[1]);
1546 
1547  return retVal;
1548 }
1549 
1550 //---------------------------------------------------------
1551 TH2D** StEStructSupport::buildChargeTypesSumOfRatios(const char* name, int opt, float *sf){
1552  TH2D** retVal= buildChargeTypesSumOfRatios(name, opt, sf, 0);
1553  if (!retVal) {
1554  return retVal;
1555  }
1556 
1557  float *nPairs = getChargeNumber(0);
1558  for (int iType=0;iType<4;iType++) {
1559  retVal[iType]->Scale(nPairs[iType]);
1560 
1561  }
1562  for (int iType=0;iType<4;iType++) {
1563  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
1564  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
1565  }
1566  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
1567  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
1568  }
1569  }
1570  for (int iz=1;iz<mNumZBins;iz++) {
1571  TH2D** tmpVal= buildChargeTypesSumOfRatios(name, opt, sf, iz);
1572  float *tempPairs = getChargeNumber(iz);
1573  for (int iType=0;iType<4;iType++) {
1574  retVal[iType]->Add(tmpVal[iType],tempPairs[iType]);
1575  nPairs[iType]+=tempPairs[iType];
1576  delete tmpVal[iType];
1577  }
1578  delete [] tmpVal;
1579  delete [] tempPairs;
1580  }
1581  for (int iType=0;iType<4;iType++) {
1582  if (nPairs[iType]) {
1583  retVal[iType]->Scale(1.0/nPairs[iType]);
1584  } else {
1585  retVal[iType]->Scale(0.0);
1586  }
1587  }
1588  delete [] nPairs;
1589  return retVal;
1590 }
1591 //---------------------------------------------------------
1592 TH2D** StEStructSupport::buildChargeTypesSumOfRatios(const char* name, int opt, float* sf, int zBin){
1593  // finds LS and US same as buildChargeTypes, but here CI and CD are sums of ratios
1594 
1595  // build hist types = LS, US, CD, CI
1596  // eight input histograms ++,+-,-+,-- for Sib and Mix
1597  TH2D** hlocal=getLocalClones(name,zBin);
1598  if (!hlocal) {
1599  return hlocal;
1600  }
1601 
1602  if(mnpairs){ // manual scaling
1603  for(int i=0;i<8;i++) hlocal[i]->Scale(1.0/mnpairs[i]);
1604  } else rescale(hlocal,zBin); // automatic scaling, norm sib to pairs per event, norm mix to sib
1605 
1606  if(strstr(name,"DEta") || strstr(name,"SEta"))fixDEtaGeometric((TH2**)hlocal,8); // does nothing unless mapplyDEtaFix is set
1607 
1608  // four returned hists
1609  TH2D** retVal= new TH2D*[4]; // 0=LS 1=US 2=CD=LS-US 3=CI=LS+US
1610  const char* oldName[4]={"Sibpp","Sibpm","Sibmp","Sibmm"};
1611  const char* newName[4]={"LS ","US ","CD ","CI "};
1612  const char* oldTitle[4]={"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-"};
1613  const char* newTitle[4]={"LS : ++ + -- ","US : +- + -+ ","CD: ++ + -- - +- - -+ ","CI : ++ + -- + +- + -+ "};
1614  for(int i=0;i<4;i++){
1615  retVal[i]=(TH2D*)hlocal[i]->Clone();
1616  retVal[i]->SetName(swapIn(hlocal[i]->GetName(),oldName[i],newName[i]));
1617  retVal[i]->SetTitle(swapIn(hlocal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1618  retVal[i]->Scale(0.); // zero the hists
1619  }
1620 
1621  hlocal[0]->Add(hlocal[3]); // ++ + -- Sib
1622  if(hlocal[2]) hlocal[1]->Add(hlocal[2]); // +- + -+ Sib, backward compatible
1623  hlocal[4]->Add(hlocal[7]); // ++ + -- Mix
1624  if(hlocal[6]) hlocal[5]->Add(hlocal[6]); // +- + -+ Mix
1625 
1626  float scf1=0.;
1627  float scf2=0.;
1628  if(sf){
1629  scf1=sf[0];
1630  scf2=sf[1];
1631  }
1632  // if requested, scale bg to require correlations>=0 where statistics are large
1633  if(1==mbgMode){
1634  scaleBackGround(hlocal[0],hlocal[4],scf1);
1635  scaleBackGround(hlocal[1],hlocal[5],scf2);
1636  if(opt==3)scaleBackGround(hlocal[3],hlocal[7]);
1637  }
1638 
1639  retVal[0]->Add(hlocal[0]); // LS sib
1640  retVal[1]->Add(hlocal[1]); // US sib
1641 
1642  // opt: 0 = delta-rho/rho_mix; 1 = delta-rho; 2 = delta-rho/sqrt(rho_mix)
1643 
1644  for(int i=0;i<2;i++){ // loop over LS and US
1645  retVal[i]->Add(hlocal[i+4],-1.0); // delta-rho
1646  }
1647  retVal[2]->Add(retVal[0],retVal[1],1.,-1.); // CD sib
1648  retVal[3]->Add(retVal[0],retVal[1],1.,1.); // CI sib
1649 
1650  if(opt!=1) { // retVal has LS, US, CD, CI sib; use hlocal to take ratios
1651  hlocal[6]->Add(hlocal[4], hlocal[5]);
1652  hlocal[7]->Add(hlocal[4], hlocal[5]); // hlocal[4] starts mix for denominators
1653 
1654  for(int i=0;i<4;i++) { //2
1655  retVal[i]->Divide(hlocal[i+4]); // delta-rho/rho_mix
1656  // NOTE: now CI = (LS_sib + US_sib) / (LS_mix + US_mix), not sum of ratios
1657 
1658  if(opt>=2){ // dN/d\eta * (delta-rho/rho_mix)
1659  TString hSibName("NEventsSib_zBuf_"); hSibName += zBin;
1660  TH1 *hNum; mtf->GetObject(hSibName.Data(),hNum);
1661  double dNdeta = hNum->GetMean()/2;
1662  retVal[i]->Scale(dNdeta/(2*3.1415926));
1663  }
1664  }
1665  }
1666 
1667  // Free local histograms.
1668  for(int i=0;i<_pair_totalmax;i++) {
1669  delete hlocal[i];
1670  }
1671  delete hlocal;
1672 
1673  return retVal;
1674 }
1675 //---------------------------------------------------------
1676 // When we have imposed some type of pt cuts we may have introduced
1677 // a covariance which will show up here. Subtract mixed pairs
1678 // to reduce this contribution.
1679 TH2D** StEStructSupport::buildPtChargeTypes(const char* name, int opt, int subtract){
1680  // See comments in buildChargeTypes.
1681  TH2D** retVal= buildPtCommon(name, opt, subtract);
1682  if (!retVal) {
1683  return retVal;
1684  }
1685 
1686  // Fix up name and title of histograms we return
1687  for (int iType=0;iType<4;iType++) {
1688  if (strstr(retVal[iType]->GetName(),"_zBuf_0")) {
1689  retVal[iType]->SetName(swapIn(retVal[iType]->GetName(),"_zBuf_0",""));
1690  }
1691  if (strstr(retVal[iType]->GetTitle(),"_zBuf_0")) {
1692  retVal[iType]->SetTitle(swapIn(retVal[iType]->GetTitle(),"_zBuf_0",""));
1693  }
1694  }
1695  const char* oldName[4] = {"PP ","PM ","MP ","MM "};
1696  const char* newName[4] = {"LS ","US ","CD ","CI "};
1697  const char* oldTitle[4] = {"PP : ++ ","PM : +- ","MP: -+ ","MM : -- "};
1698  const char* newTitle[4] = {"LS : ++ + -- ","US : +- + -+ ","CD: ++ + -- - +- - -+ ","CI : ++ + -- + +- + -+ "};
1699  for (int i=0;i<4;i++) {
1700  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
1701  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1702  }
1703 
1704  // Form charge combinations LS, US, CD and CI.
1705  // Need to scale PP by (n+/Nch)^2, PM by n+n-/Nch^2 and MM by (n-/Nch)^2 to combine terms properly
1706  double *scale = getScaleFactors();
1707  for (int iType=0;iType<4;iType++) {
1708  // Need to undo scaling with charge specific normalization and redo with CI
1709  if (0 == opt) {
1710  retVal[iType]->Scale(scale[iType]/pow(scale[4],2));
1711  } else if (1 == opt) {
1712  retVal[iType]->Scale(1.0);
1713  } else if (2 == opt) {
1714  retVal[iType]->Scale(sqrt(scale[iType])/scale[4]);
1715  }
1716  }
1717  delete [] scale;
1718  retVal[0]->Add(retVal[3]);
1719  if (mIdenticalPair) {
1720  retVal[1]->Scale(2);
1721  } else {
1722  retVal[1]->Add(retVal[2]);
1723  }
1724  retVal[2]->Add(retVal[0],retVal[1],1,-1);
1725  retVal[3]->Add(retVal[0],retVal[1]);
1726  return retVal;
1727 }
1728 //---------------------------------------------------------
1729 // This routine is not used except for checking zBuffer dependence.
1730 TH2D** StEStructSupport::buildPtChargeTypes(const char* name, int opt, int subtract, int zBin) {
1731 
1732  TH2D** retVal= buildPtCommon(name, opt, subtract, zBin);
1733  if (!retVal) {
1734  return retVal;
1735  }
1736 
1737  // Fix up name and title of histograms we return
1738  const char* oldName[4] = {"Sibpp","Sibpm","Sibmp","Sibmm"};
1739  const char* newName[4] = {"LS ","US ","CD ","CI "};
1740  const char* oldTitle[4] = {"Sibling : +.+","Sibling : +.-","Sibling : -.+","Sibling : -.-"};
1741  const char* newTitle[4] = {"LS : ++ + -- ","US : +- + -+ ","CD: ++ + -- - +- - -+ ","CI : ++ + -- + +- + -+ "};
1742  for(int i=0;i<4;i++){
1743  retVal[i]->SetName(swapIn(retVal[i]->GetName(),oldName[i],newName[i]));
1744  retVal[i]->SetTitle(swapIn(retVal[i]->GetTitle(),oldTitle[i],newTitle[i]));
1745  }
1746 
1747  // Form charge combinations.
1748  retVal[0]->Add(retVal[3]);
1749  retVal[1]->Add(retVal[2]);
1750  retVal[2]->Add(retVal[0],retVal[1],1,-1);
1751  retVal[3]->Add(retVal[0],retVal[1]);
1752  return retVal;
1753 }
1754 
1755 //---------------------------------------------------------
1756 // djp I don't understand the logic here. If any mixed bin exceeds the sibling
1757 // bin by three sigma and is at least 45% of the maximum sibling bin then
1758 // we scale entire mixed histogram by sibling/mixed for this one bin.
1759 // And if it happens for more than one bin we use the last one!
1760 //
1761 // Finally noticed, when we explicitly give a scale factor we use that.
1762 void StEStructSupport::scaleBackGround(TH2D* sib, TH2D* mix, float sf) {
1763  float alpha=1.0;
1764  if (sf!=0.) {
1765  alpha=sf;
1766  } else {
1767  float hmax = mix->GetMaximum();
1768  for (int ix=1;ix<=mix->GetNbinsX();ix++) {
1769  for(int iy=1;iy<=mix->GetNbinsY();iy++) {
1770  float mval = mix->GetBinContent(ix,iy);
1771  if (mval>0.) {
1772  float sval = sib->GetBinContent(ix,iy);
1773  if (sval==0.) {
1774  continue;
1775  }
1776  float eval = sib->GetBinError(ix,iy)/sval;
1777  float dval = fabs(sval-mval)/sval;
1778  if (sval<mval && dval>3.0*eval && mval/hmax>0.45) {
1779 // if (sval/hmax>0.5) {
1780  float rval = sval/mval;
1781  if (rval<alpha) {
1782  alpha = rval;
1783  }
1784  }
1785  }
1786  }
1787  }
1788  }
1789  mix->Scale(alpha);
1790  if (!msilent) {
1791  cout<<"Scaling "<<mix->GetName()<<" by "<<alpha<<endl;
1792  }
1793 }
1794 
1795 //----------------------------------------------------------------
1796 void StEStructSupport::fixDEtaGeometric(TH2** h, int numHists) {
1797 
1798  // That other version of fixDEta is actually trying to do a complicated
1799  // acceptance correction. The geometric part is well defined, and we
1800  // do only that here.
1801  // No need to assume max and min eta are same distance from 0.
1802  // Makes it look slightly more complicated.
1803  // Take geometry from first histogram.
1804 
1805  if(!mapplyDEtaFix) return;
1806 
1807  TAxis *x = h[0]->GetXaxis();
1808  TAxis *y = h[0]->GetYaxis();
1809  double etaMax = x->GetXmax();
1810  double etaMin = x->GetXmin();
1811  double H = 0.5*( etaMax - etaMin );
1812  double x0 = 0.5*( etaMax + etaMin );
1813  for (int ix=1;ix<=x->GetNbins();ix++) {
1814  double eta1 = x->GetBinLowEdge(ix);
1815  double eta2 = x->GetBinUpEdge(ix);
1816  double corr;
1817  if (eta2 <= x0) {
1818  corr = H / ( 0.5*(eta2+eta1) - etaMin );
1819  } else if (eta1 >= x0) {
1820  corr = H / ( etaMax - 0.5*(eta2+eta1) );
1821  } else {
1822  double area = H*(eta2-eta1);
1823  double aplus = (eta2-x0)*0.5*(H + etaMax-eta2);
1824  double aminus = (x0-eta1)*0.5*(H + eta1-etaMin);
1825  corr = area / (aplus + aminus);
1826  }
1827  for (int iy=1;iy<=y->GetNbins();iy++) {
1828  for (int ih=0;ih<numHists;ih++) {
1829  double val = h[ih]->GetBinContent(ix,iy);
1830  double err = h[ih]->GetBinError(ix,iy);
1831  h[ih]->SetBinContent(ix,iy,val*corr);
1832  h[ih]->SetBinError(ix,iy,err*corr);
1833  }
1834  }
1835  }
1836 }
1837 
1838 //----------------------------------------------------------------
1839 void StEStructSupport::fixDEta(TH2** h, int numHists) {
1840 
1841  // Expect histograms in groups of eight. The first eight are number,
1842  // with the first four being sibling and the second four being
1843  // mixed. If there are more than eight they will be pt weighted.
1844  // I find acceptance in the pair hole region is different for ++ and
1845  // +-, so we create two acceptance histograms, ++ + -- and +- + -+
1846  // using the mixed number. Not sure how to normalize the acceptance
1847  // histograms. For now try setting acceptance for \eta_\Delta = 0
1848  // to 1, excluding the pair cut hole.
1849 
1850  if(!mapplyDEtaFix) return;
1851 
1852  TH2D *acc[2];
1853  acc[0] = (TH2D *) h[4]->Clone();
1854  acc[0]->Add(h[7]);
1855  acc[1] = (TH2D *) h[5]->Clone();
1856  acc[1]->Add(h[6]);
1857  int n = acc[0]->GetNbinsX();
1858  double scale[] = {0, 0};
1859  if (n%2) {
1860  int mid = 1 + n/2;
1861  for (int iy=12;iy<=25;iy++) {
1862  scale[0] += acc[0]->GetBinContent(mid,iy);
1863  scale[1] += acc[1]->GetBinContent(mid,iy);
1864  }
1865  scale[0] /= 14;
1866  scale[1] /= 14;
1867  } else {
1868  int mid = n/2;
1869  for (int iy=12;iy<=25;iy++) {
1870  scale[0] += acc[0]->GetBinContent(mid,iy);
1871  scale[0] += acc[0]->GetBinContent(mid+1,iy);
1872  scale[1] += acc[1]->GetBinContent(mid,iy);
1873  scale[1] += acc[1]->GetBinContent(mid+1,iy);
1874  }
1875  scale[0] /= 28;
1876  scale[1] /= 28;
1877  }
1878  if (scale[0] > 0) {
1879  acc[0]->Scale(1.0/scale[0]);
1880  }
1881  if (scale[1] > 0) {
1882  acc[1]->Scale(1.0/scale[1]);
1883  }
1884 
1885  for (int ig=0;ig<numHists;ig+=8) {
1886  if (acc[0]->Integral() > 0) {
1887  h[ig+0]->Divide(acc[0]);
1888  h[ig+1]->Divide(acc[1]);
1889  h[ig+2]->Divide(acc[1]);
1890  h[ig+3]->Divide(acc[0]);
1891  h[ig+4]->Divide(acc[0]);
1892  h[ig+5]->Divide(acc[1]);
1893  h[ig+6]->Divide(acc[1]);
1894  h[ig+7]->Divide(acc[0]);
1895  }
1896  }
1897  delete acc[0];
1898  delete acc[1];
1899 }
1900 
1901 //---------------------------------------------------------
1902 void StEStructSupport::writeAscii(TH2D** h, int iHist, const char* fname, int optErrors){
1903 
1904  ofstream of(fname);
1905  int xbins=h[iHist]->GetNbinsX();
1906  TAxis * xa= h[iHist]->GetXaxis();
1907  of<<"# Histogram Name = "<<h[iHist]->GetName()<<endl;
1908  of<<"# X-axis: min="<<xa->GetBinLowEdge(1)<<" max="<<xa->GetBinUpEdge(xbins)<<" nbins="<<xbins<<endl;
1909 
1910  int ybins=h[iHist]->GetNbinsY();
1911  TAxis * ya= h[iHist]->GetYaxis();
1912  if(ybins>1){
1913  of<<"# Y-axis: min="<<ya->GetBinLowEdge(1)<<" max="<<ya->GetBinUpEdge(ybins)<<" nbins="<<ybins<<endl;
1914  of<<"# ix iy Value "<<endl;
1915  for(int i=1;i<=xbins;i++){
1916  for(int j=1;j<=ybins;j++){
1917  of<<i<<" "<<j<<" "<<h[iHist]->GetBinContent(i,j);
1918  if(optErrors>0)of<<" "<<h[iHist]->GetBinError(i,j);
1919  of<<endl;
1920  }
1921  }
1922  of<<endl;
1923  } else {
1924  of<<"# ix Value "<<endl;
1925  for(int i=1;i<=xbins;i++){
1926  of<<i<<" "<<h[iHist]->GetBinContent(i);
1927  if(optErrors>0)of<<" "<<h[iHist]->GetBinError(i);
1928  of<<endl;
1929  }
1930  of<<endl;
1931  }
1932 
1933  of.close();
1934 }
1935 
1936 //---------------------------------------------------------
1937 char* StEStructSupport::swapIn(const char* name, const char* s1, const char* s2){
1938 
1939  // looks for s1 in name and replaces with s2
1940  //
1941  // in perl $name=~s/$s1/$s2/;
1942 
1943  if(mtmpString) delete [] mtmpString;
1944  mtmpString=NULL;
1945 
1946  if(!name) return mtmpString;
1947 
1948  int len=strlen(name);
1949  char* tmp=new char[len+1];
1950  strcpy(tmp,name);
1951 
1952  char* ptr1;
1953  if(!(ptr1=strstr(tmp,s1))) {
1954  delete [] tmp;
1955  return mtmpString;
1956  }
1957 
1958  int len1=strlen(s1);
1959  mtmpString=new char[256];
1960 
1961  std::ostringstream ts;
1962  char* ptr2=ptr1;
1963  *ptr1='\0';
1964  ts<<tmp;
1965  if(s2)ts<<s2;
1966  ptr1+=len1;
1967  ts<<ptr1;
1968  *ptr2=' ';
1969 
1970  delete [] tmp;
1971  strcpy(mtmpString,(ts.str()).c_str());
1972  return mtmpString;
1973 }
1974 
1975 /***********************************************************************
1976  *
1977  * $Log: StEStructSupport.cxx,v $
1978  * Revision 1.30 2012/11/16 21:27:23 prindle
1979  * AS, SS histograms. Get eta limits from histogram in file. Check for empty bins in ratio (showed up as offset in correlations). Scale histograms according to actual number of pairs, don't assume 1/2 for example.
1980  *
1981  * Revision 1.29 2012/06/11 14:35:33 fisyak
1982  * std namespace
1983  *
1984  * Revision 1.28 2011/08/02 20:42:24 prindle
1985  * Added YtYtVolumeNormalization.
1986  * Fixed calculation of error for YtYt \Delta\rho/sqrt(\rho_{ref})
1987  * Added error calculation for p_t histograms
1988  * Added warning when either \rho_{sib} or \rho_{ref} has an empty bin. Set ratio
1989  * bin to 0 instead of -1.
1990  *
1991  * Revision 1.27 2010/09/02 21:31:14 prindle
1992  * Support: Can't find evidence that YtYt correlation depends on z vertex.
1993  * Add numerators and denominators then take ratio. Need a new
1994  * rescale method independent of z bin. Looks like we can normalize
1995  * mixed so \int{sib/mix} = number of bins (what we have recently been
1996  * doing) or \int{sib} = \int{mix} and the former is more snesitive
1997  * to bins with very few counts. That isn't important for angular
1998  * histograms but is for (yt,yt). I am calling the \int{sib} = \int(mix}
1999  * Yt normalization (even though it is what we did long ago).
2000  *
2001  * Revision 1.26 2010/06/23 22:33:50 prindle
2002  * In HAdd we distinguish between the parent distributions of the
2003  * two particles.
2004  * In Support I fixed a number of problems in the Pt correlation section.
2005  *
2006  * Revision 1.25 2010/03/02 21:48:30 prindle
2007  * Fix addDensities (for checking pair cuts)
2008  * Lots of small changes
2009  *
2010  * Revision 1.24 2009/11/09 21:32:59 prindle
2011  * Fix warnings about casting char * to a const char * by redeclaring as const char *.
2012  *
2013  * Revision 1.23 2009/08/17 23:05:20 dkettler
2014  * Normalization fix
2015  *
2016  * Revision 1.21 2009/05/08 00:21:42 prindle
2017  * In StEStructHadd remove support for old style of histogram names, do a better job calculating
2018  * errors (at least for number (\eta_\Delta,\phi_\Delta) histograms), double bins which
2019  * have an edge in the center (purely cosmetic when looking at intermediate histograms).
2020  * In StEStructSupport check for existance of histograms and return gracefully.
2021  * Code in buildChargeTypes and buildPtChargeTypes was essentially duplicate of code
2022  * in buildCommon and buildPtCommon so I refactored to reduce redundancy.
2023  *
2024  * Revision 1.20 2009/02/03 14:30:23 fisyak
2025  * Add missing includes for ROOT 5.22
2026  *
2027  * Revision 1.19 2008/12/02 23:52:53 prindle
2028  * Get information about histogram XX being symmetrized from CutBin.
2029  * Changed TH1* to TH2D* in many places hoping to be able to plot DEtaDPhi
2030  * as colz (doesn't work yet).
2031  * Added direct calculation of \Delta\rho/\rho_{ref} (and similar) which is
2032  * needed for YtYt correlations.
2033  *
2034  * Revision 1.18 2008/03/19 22:08:38 prindle
2035  * Use GetObject instead of Get for type safety. Stop deleting objects we didn't create.
2036  * Treat \Delta\rho = d^2n/dEtadphi (rho - rho_ref)/rho_ref as basic unit when combining
2037  * centralities and z bins.
2038  *
2039  * This code should be check in more detail before being completely relied upon.
2040  *
2041  * Revision 1.17 2007/11/26 20:07:19 prindle
2042  * Modified to average \Delta\rho/sqrt(\rho) over z-bins (if more than one z-bin
2043  * present for given centrality. Note: I weight by number of tracks, not number of
2044  * pairs. This is important when we are also combining different centralities (which
2045  * I do by combining centrality tag with z-bin tag in macro/addCentralities.)
2046  *
2047  * Scale mixed histograms by number of events. Integral of \Delta\rho need not be 0.
2048  *
2049  * delete items that are created and valgrind complained were lost. (Not a big deal
2050  * since macro is run once)
2051  *
2052  * [Still need to commit StEStructHAdd.cxx which cvs complained that check-update failed.]
2053  *
2054  * Revision 1.16 2007/05/27 22:46:01 msd
2055  * Added buildChargeTypes mode 3 which takes rho_ref from track counts.
2056  * Added buildChargeTypeSumOfRatios.
2057  *
2058  * Revision 1.15 2007/01/26 17:20:58 msd
2059  * Updated HAdd for new binning scheme.
2060  * Improved Support::buildChargeTypes.
2061  *
2062  * Revision 1.14 2006/12/14 20:07:11 prindle
2063  * I was calculating \Delta\rho/sqrt(rho) for ++, +-, -+ and --
2064  * and then combining those into LS, US, CD and CI. The was wrong
2065  * and now I am doing it correctly. For CI this makes only a slight
2066  * change, it seems the amplitude is decreased a little. For CD
2067  * this is a bigger change. I left the old versions (with _Old appended)
2068  * for now.
2069  *
2070  * Revision 1.13 2006/10/27 00:05:32 prindle
2071  * Modified buildChargeTypes to handle case where the -+ histogram is
2072  * empty. Also tried making buildCommonTypes work, but there one of the
2073  * output histograms was intended to be -+ and most of the time that
2074  * will be empty.
2075  *
2076  * Revision 1.11 2006/04/26 18:52:12 dkettler
2077  *
2078  * Added reaction plane determination for the analysis
2079  *
2080  * Added reaction plane angle calculation
2081  *
2082  * Case 3 in buildPtChargeTypes needs to be corrected
2083  *
2084  * Revision 1.10 2006/04/25 21:04:39 msd
2085  * Added Jeff's patch for getHists to create doubles instead of floats
2086  *
2087  * Revision 1.9 2006/04/06 01:09:49 prindle
2088  * Calculating pt for each cut bin caused changes in HAdd.
2089  * The splitting of +- into +- and -+ caused changes in Support.
2090  *
2091  * Revision 1.8 2006/04/04 22:14:10 porter
2092  * fixdeta is now NOT default but included in StEStruct2ptCorrelations
2093  *
2094  * Revision 1.7 2005/09/14 17:25:37 msd
2095  * minor tweak
2096  *
2097  * Revision 1.6 2005/09/07 20:26:16 prindle
2098  *
2099  *
2100  * Support: Fixed some meory leaks.
2101  *
2102  * Revision 1.4 2005/03/08 21:56:42 porter
2103  * fixed bug in StEStructHAdd.cxx and added diagnostic option in ptcorrelations to
2104  * view individual terms separately
2105  *
2106  * Revision 1.3 2005/03/08 20:16:34 msd
2107  * included <sstream> library
2108  *
2109  * Revision 1.2 2005/03/03 01:33:05 porter
2110  * Added pt-correlations method to support and included
2111  * these histograms to the HAdd routine
2112  *
2113  * Revision 1.1 2004/07/01 00:37:17 porter
2114  * new code previously my StEStructHelper. Takes hists from correltation
2115  * pass and builds final ressults. Also the StEStructHAdd.h is a simple
2116  * replacemnt for my sumyt.C macro which could be expanded later as needed.
2117  *
2118  *
2119  *
2120  *********************************************************************/
2121 
2122 
StEStructSupport::mbgMode
int mbgMode
for normalization comparing different cuts
Definition: StEStructSupport.h:31
Generated by   doxygen 1.8.5