StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructSigmas.cxx
1 
2 #include "StEStructSigmas.h"
3 
4 #include "TMath.h"
5 #include "TH1F.h"
6 #include "TH2F.h"
7 #include "TFile.h"
8 
9 
10 ClassImp(StEStructSigmas)
11 
12 //--------------------------------------------------------------------------
13 StEStructSigmas::StEStructSigmas( char *key,
14  int nPhi, int nEta,
15  float EtaMin, float EtaMax,
16  const char *preFix ) {
17  mKey = strdup(key);
18  mNPhiBins = nPhi;
19  mNEtaBins = nEta;
20  mEtaMin = EtaMin;
21  mEtaMax = EtaMax;
22  mpreFix = strdup(preFix);
23  mEtaSumMode = 1;
24  mPhiSumMode = 1;
25  initHistograms();
26 }
27 //--------------------------------------------------------------------------
28 StEStructSigmas::~StEStructSigmas() {
29  deleteHistograms();
30  free(mKey);
31  free(mpreFix);
32 }
33 
34 void StEStructSigmas::fillHistograms() {
35  NHistograms();
36  PHistograms();
37  PNHistograms();
38 }
39 
40 
41 //--------------------------------------------------------------------------
42 void StEStructSigmas::NHistograms() {
43  char buffer[255];
44  binTupleStruct bTuple;
45  scaleTupleStruct sTuple;
46  sumTupleStruct STuple;
47 
48  TH2F *hnBins = (TH2F *) gDirectory->Get("nBins");
49  TH2F *hoffset = (TH2F *) gDirectory->Get("offset");
50  TH2F *hfUnique = (TH2F *) gDirectory->Get("fUnique");
51 
52  // We have accumulated over events.
53  // For every bin at each scale we have summed up all the small bins.
54  // Now accumulate sum, n, n^2 for all bins.
55 
56  NSig->Reset();
57  NDel->Reset();
58  NPlus->Reset();
59  NMinus->Reset();
60  NPlusMinus->Reset();
61 
62  NSigCorrection->Reset();
63  NDelCorrection->Reset();
64  NPlusCorrection->Reset();
65  NMinusCorrection->Reset();
66  NPlusMinusCorrection->Reset();
67 
68  NSigErrors->Reset();
69  NDelErrors->Reset();
70  NPlusErrors->Reset();
71  NMinusErrors->Reset();
72  NPlusMinusErrors->Reset();
73 
74  TH1D *hTotEvents[5];
75  TH1D *hNSum[16];
76  TH1D *hNDiff[2];
77  TH1D *hNPlus[5];
78  TH1D *hNMinus[5];
79  TH1D *hNPlusMinus[8];
80 
81  sprintf(buffer,"%sTotalEvents_%s", mpreFix, mKey);
82  hTotEvents[0] = (TH1D *) gDirectory->Get(buffer);
83  sprintf(buffer,"%sTotalSumEvents_%s", mpreFix, mKey);
84  hTotEvents[1] = (TH1D *) gDirectory->Get(buffer);
85  sprintf(buffer,"%sTotalPlusEvents_%s", mpreFix, mKey);
86  hTotEvents[2] = (TH1D *) gDirectory->Get(buffer);
87  sprintf(buffer,"%sTotalMinusEvents_%s", mpreFix, mKey);
88  hTotEvents[3] = (TH1D *) gDirectory->Get(buffer);
89  sprintf(buffer,"%sTotalPlusMinusEvents_%s", mpreFix, mKey);
90  hTotEvents[4] = (TH1D *) gDirectory->Get(buffer);
91 
92  for (int jStat=0;jStat<16;jStat++) {
93  sprintf( buffer,"%sNSum%s_%i", mpreFix, mKey, jStat );
94  hNSum[jStat] = (TH1D *) gDirectory->Get(buffer);
95  }
96  for (int jStat=0;jStat<2;jStat++) {
97  sprintf( buffer,"%sNDiff%s_%i", mpreFix, mKey, jStat );
98  hNDiff[jStat] = (TH1D *) gDirectory->Get(buffer);
99  }
100  for (int jStat=0;jStat<5;jStat++) {
101  sprintf( buffer,"%sNPlus%s_%i", mpreFix, mKey, jStat );
102  hNPlus[jStat] = (TH1D *) gDirectory->Get(buffer);
103  }
104  for (int jStat=0;jStat<5;jStat++) {
105  sprintf( buffer,"%sNMinus%s_%i", mpreFix, mKey, jStat );
106  hNMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
107  }
108  for (int jStat=0;jStat<8;jStat++) {
109  sprintf( buffer,"%sNPlusMinus%s_%i", mpreFix, mKey, jStat );
110  hNPlusMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
111  }
112 
113  double vecCI[NPHIBINS][NETABINS][3];
114  double vecCD[NPHIBINS][NETABINS][3];
115  double vecPP[NPHIBINS][NETABINS][3];
116  double vecMM[NPHIBINS][NETABINS][3];
117  double vecPM[NPHIBINS][NETABINS][3];
118  double BCI = 0, nBCI = 0, BCD = 0, nBCD = 0, BPP = 0, nBPP = 0;
119  double BMM = 0, nBMM = 0, BPM = 0, nBPM = 0;
120  int mBin, oBin;
121  for (int iPhi=1;iPhi<=mNPhiBins;iPhi++) {
122  int jPhi = iPhi - 1;
123  for (int iEta=1;iEta<=mNEtaBins;iEta++) {
124  int jEta = iEta - 1;
125  for (int i=0;i<3;i++) {
126  vecCI[jPhi][jEta][i] = 0;
127  vecCD[jPhi][jEta][i] = 0;
128  vecPP[jPhi][jEta][i] = 0;
129  vecMM[jPhi][jEta][i] = 0;
130  vecPM[jPhi][jEta][i] = 0;
131  }
132 
133  double NS[16], ND[2], NP[5], NM[5], NC[8];
134  double sigSq, psigSq, msigSq;
135  double Ss = 0, nSs = 0, DSs = 0, nDSs = 0;
136  double Sd = 0, nSd = 0, DSd = 0, nDSd = 0;
137  double Sp = 0, nSp = 0, DSp = 0, nDSp = 0;
138  double Sm = 0, nSm = 0, DSm = 0, nDSm = 0;
139  double Sc = 0, nSc = 0, DSc = 0, nDSc = 0;
140  double sumEvents, plusEvents, minusEvents, plusMinusEvents;
141 
142  mBin = (int) hnBins->GetBinContent(iPhi,iEta);
143  oBin = (int) hoffset->GetBinContent(iPhi,iEta);
144 
145  int iEtaBin = 0, iPhiBin = 1;
146  double mCI[2][2], mCD[2][2], mPP[2][2], mMM[2][2], mPM[2][2];
147  double vCI[2], vCD[2], vPP[2], vMM[2], vPM[2];
148  double mInv[2][2], mDet;
149  int nEtaBin = getNumEtaBins(iEta);
150 // int nPhiBin = getNumPhiBins(iEta);
151  float phi1, phi2, eta1, eta2;
152  double f3;
153  for (int i=0;i<2;i++) {
154  vCI[i] = 0;
155  vCD[i] = 0;
156  vPP[i] = 0;
157  vMM[i] = 0;
158  vPM[i] = 0;
159  for (int j=0;j<2;j++) {
160  mCI[i][j] = 0;
161  mCD[i][j] = 0;
162  mPP[i][j] = 0;
163  mMM[i][j] = 0;
164  mPM[i][j] = 0;
165  }
166  }
167  for (int iBin=oBin+1;iBin<=oBin+mBin;iBin++) {
168  iEtaBin++;
169  if (iEtaBin > nEtaBin) {
170  iEtaBin = 1;
171  iPhiBin++;
172  }
173  phi1 = 0 + 2*3.1415926*getPhiStart(iPhiBin,iPhi)/mNPhiBins;
174  phi2 = phi1 + iPhi * 2*3.1415926 / mNPhiBins;
175  eta1 = mEtaMin + (mEtaMax-mEtaMin)*getEtaStart(iEtaBin,iEta)/mNEtaBins;
176  eta2 = eta1 + iEta * (mEtaMax-mEtaMin) / mNEtaBins;
177  f3 = (pow(eta2,3)-pow(eta1,3))/(3*(eta2-eta1));
178  sumEvents = hTotEvents[0]->GetBinContent(iBin);
179  if (sumEvents > 0) {
180  for (int jStat=0;jStat<16;jStat++) {
181  NS[jStat] = hNSum[jStat]->GetBinContent(iBin) / sumEvents;
182  }
183  for (int jStat=0;jStat<2;jStat++) {
184  ND[jStat] = hNDiff[jStat]->GetBinContent(iBin) / sumEvents;
185  }
186  double NSum = NS[0] + NS[1];
187  if (NSum > 0) {
188  double SumSq = NS[2] + 2*NS[3] + NS[4];
189  sigSq = (SumSq - NSum*NSum) / NSum;
190  Ss += sigSq - 1;
191  nSs += 1;
192  DSs += (4 + sigSq/NSum) * sigSq / sumEvents;
193  nDSs += sqrt(hfUnique->GetBinContent(iPhi,iEta));
194  bTuple.type = 1;
195  bTuple.phiScale = iPhi;
196  bTuple.etaScale = iEta;
197  bTuple.phi = (phi1+phi2)/2;
198  bTuple.eta = (eta1+eta2)/2;
199  bTuple.sig2 = sigSq-1;
200  bTuple.sig2_1 = 0;
201  bTuple.sig2_2 = 0;
202  bTuple.nbar = NSum;
203  bTuple.events = sumEvents;
204  bTuple.f3 = f3;
205  binTuple->Fill(&bTuple.type);
206  vecCI[jPhi][jEta][0] += 1;
207  vecCI[jPhi][jEta][1] += sigSq-1;
208  vecCI[jPhi][jEta][2] += f3;
209  if (nEtaBin > 3) {
210  mCI[0][0] += 1;
211  mCI[0][1] += f3;
212  mCI[1][0] += f3;
213  mCI[1][1] += f3*f3;
214  vCI[0] += sigSq-1;
215  vCI[1] += (sigSq-1)*f3;
216  }
217 
218  double NDiff = NS[0] - NS[1];
219  double DiffSq = NS[2] - 2*NS[3] + NS[4];
220  sigSq = (DiffSq - NDiff*NDiff) / NSum;
221  Sd += sigSq - 1;
222  nSd += 1;
223  DSd += (4 + sigSq/NSum) * sigSq / sumEvents;
224  nDSd += sqrt(hfUnique->GetBinContent(iPhi,iEta));
225  bTuple.type = 2;
226  bTuple.phiScale = iPhi;
227  bTuple.etaScale = iEta;
228  bTuple.phi = (phi1+phi2)/2;
229  bTuple.eta = (eta1+eta2)/2;
230  bTuple.sig2 = sigSq-1;
231  bTuple.sig2_1 = 0;
232  bTuple.sig2_2 = 0;
233  bTuple.nbar = NSum;
234  bTuple.events = sumEvents;
235  bTuple.f3 = f3;
236  binTuple->Fill(&bTuple.type);
237  vecCD[jPhi][jEta][0] += 1;
238  vecCD[jPhi][jEta][1] += sigSq-1;
239  vecCD[jPhi][jEta][2] += f3;
240  if (nEtaBin > 3) {
241  mCD[0][0] += 1;
242  mCD[0][1] += f3;
243  mCD[1][0] += f3;
244  mCD[1][1] += f3*f3;
245  vCD[0] += sigSq-1;
246  vCD[1] += (sigSq-1)*f3;
247  }
248  }
249  }
250 
251  psigSq = 0;
252  plusEvents = hTotEvents[0]->GetBinContent(iBin);
253  if (plusEvents > 0) {
254  for (int jStat=0;jStat<5;jStat++) {
255  NP[jStat] = hNPlus[jStat]->GetBinContent(iBin) / plusEvents;
256  }
257  if (NP[0] > 0) {
258  sigSq = (NP[1] - NP[0]*NP[0]) / NP[0];
259  Sp += sigSq - 1;
260  nSp += 1;
261  DSp += (4 + sigSq/NP[0]) * sigSq / plusEvents;
262  nDSp += sqrt(hfUnique->GetBinContent(iPhi,iEta));
263  psigSq = sigSq;
264  bTuple.type = 3;
265  bTuple.phiScale = iPhi;
266  bTuple.etaScale = iEta;
267  bTuple.phi = (phi1+phi2)/2;
268  bTuple.eta = (eta1+eta2)/2;
269  bTuple.sig2 = sigSq-1;
270  bTuple.sig2_1 = 0;
271  bTuple.sig2_2 = 0;
272  bTuple.nbar = NP[0];
273  bTuple.events = plusEvents;
274  bTuple.f3 = f3;
275  binTuple->Fill(&bTuple.type);
276  vecPP[jPhi][jEta][0] += 1;
277  vecPP[jPhi][jEta][1] += sigSq-1;
278  vecPP[jPhi][jEta][2] += f3;
279  if (nEtaBin > 3) {
280  mPP[0][0] += 1;
281  mPP[0][1] += f3;
282  mPP[1][0] += f3;
283  mPP[1][1] += f3*f3;
284  vPP[0] += sigSq-1;
285  vPP[1] += (sigSq-1)*f3;
286  }
287  }
288  } else {
289  for (int jStat=0;jStat<5;jStat++) {
290  NP[jStat] = 0;
291  }
292  }
293 
294  msigSq = 0;
295  minusEvents = hTotEvents[0]->GetBinContent(iBin);
296  if (minusEvents > 0) {
297  for (int jStat=0;jStat<5;jStat++) {
298  NM[jStat] = hNMinus[jStat]->GetBinContent(iBin) / minusEvents;
299  }
300  if (NM[0] > 0) {
301  sigSq = (NM[1] - NM[0]*NM[0]) / NM[0];
302  Sm += sigSq - 1;
303  nSm += 1;
304  DSm += (4 + sigSq/NM[0]) * sigSq / minusEvents;
305  nDSm += sqrt(hfUnique->GetBinContent(iPhi,iEta));
306  msigSq = sigSq;
307  bTuple.type = 4;
308  bTuple.phiScale = iPhi;
309  bTuple.etaScale = iEta;
310  bTuple.phi = (phi1+phi2)/2;
311  bTuple.eta = (eta1+eta2)/2;
312  bTuple.sig2 = sigSq-1;
313  bTuple.sig2_1 = 0;
314  bTuple.sig2_2 = 0;
315  bTuple.nbar = NM[0];
316  bTuple.events = minusEvents;
317  bTuple.f3 = f3;
318  binTuple->Fill(&bTuple.type);
319  vecMM[jPhi][jEta][0] += 1;
320  vecMM[jPhi][jEta][1] += sigSq-1;
321  vecMM[jPhi][jEta][2] += f3;
322  if (nEtaBin > 3) {
323  mMM[0][0] += 1;
324  mMM[0][1] += f3;
325  mMM[1][0] += f3;
326  mMM[1][1] += f3*f3;
327  vMM[0] += sigSq-1;
328  vMM[1] += (sigSq-1)*f3;
329  }
330  }
331  } else {
332  for (int jStat=0;jStat<5;jStat++) {
333  NM[jStat] = 0;
334  }
335  }
336 
337 // Accumulation of hNPlusMinus required (n+ > 0) && (n- > 0)
338 // I think this is wrong. Instead if I require (n+ > 0) || (n- > 0)
339 // I can retrieve appropriate terms from NP and NM, at least for the
340 // phiPt case.
341  plusMinusEvents = hTotEvents[0]->GetBinContent(iBin);
342  if (plusMinusEvents > 0) {
343  for (int jStat=0;jStat<8;jStat++) {
344  NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
345  }
346  if (NP[0]*NM[0] > 0) {
347  sigSq = (NC[2] - NP[0]*NM[0])/sqrt(NP[0]*NM[0]);
348  Sc += sigSq;
349  nSc += 1;
350  DSc += (psigSq + msigSq) / plusMinusEvents +
351  (1/NC[0]+1/NC[1])*sigSq*sigSq/4;
352  nDSc += sqrt(hfUnique->GetBinContent(iPhi,iEta));
353  bTuple.type = 5;
354  bTuple.phiScale = iPhi;
355  bTuple.etaScale = iEta;
356  bTuple.phi = (phi1+phi2)/2;
357  bTuple.eta = (eta1+eta2)/2;
358  bTuple.sig2 = sigSq;
359  bTuple.sig2_1 = 0;
360  bTuple.sig2_2 = 0;
361  bTuple.nbar = sqrt(NP[0]*NM[0]);
362  bTuple.events = plusMinusEvents;
363  bTuple.f3 = f3;
364  binTuple->Fill(&bTuple.type);
365  vecPM[jPhi][jEta][0] += 1;
366  vecPM[jPhi][jEta][1] += sigSq;
367  vecPM[jPhi][jEta][2] += f3;
368  if (nEtaBin > 3) {
369  mPM[0][0] += 1;
370  mPM[0][1] += f3;
371  mPM[1][0] += f3;
372  mPM[1][1] += f3*f3;
373  vPM[0] += sigSq;
374  vPM[1] += sigSq*f3;
375  }
376  }
377  }
378 // plusMinusEvents = hTotEvents[0]->GetBinContent(iBin);
379 // if (plusMinusEvents > 0) {
380 // for (int jStat=0;jStat<8;jStat++) {
381 // NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
382 // }
383 // if (NC[0]*NC[1] > 0) {
384 // sigSq = (NC[2] - NC[0]*NC[1])/sqrt(NC[0]*NC[1]);
385 // Sc += sigSq;
386 // nSc += 1;
387 // DSc += (psigSq + msigSq) / plusMinusEvents +
388 // (1/NC[0]+1/NC[1])*sigSq*sigSq/4;
389 // nDSc += sqrt(hfUnique->GetBinContent(iPhi,iEta));
390 // }
391 // }
392  }
393 
394  sTuple.phiScale = iPhi;
395  sTuple.etaScale = iEta;
396  mDet = mCI[0][0]*mCI[1][1] - mCI[0][1]*mCI[1][0];
397  if (mDet > 0.00000001) {
398  mInv[0][0] = +mCI[1][1] / mDet;
399  mInv[1][0] = -mCI[1][0] / mDet;
400  mInv[0][1] = -mCI[0][1] / mDet;
401  mInv[1][1] = +mCI[0][0] / mDet;
402  sTuple.type = 1;
403  sTuple.A = (vCI[0]*mInv[0][0] + vCI[1]*mInv[1][0]) / (iPhi*iEta);
404  sTuple.B = (vCI[0]*mInv[0][1] + vCI[1]*mInv[1][1]) / (iPhi*iEta);
405  sTuple.nBins = mCI[0][0];
406  sTuple.f3 = mCI[0][1];
407  sTuple.f3sq = mCI[1][1];
408  sTuple.sig2 = vCI[0];
409  sTuple.sig2f3 = vCI[1];
410  scaleTuple->Fill(&sTuple.type);
411  BCI += sTuple.B;
412  nBCI += 1;
413  }
414  mDet = mCD[0][0]*mCD[1][1] - mCD[0][1]*mCD[1][0];
415  if (mDet > 0.00000001) {
416  mInv[0][0] = +mCD[1][1] / mDet;
417  mInv[1][0] = -mCD[1][0] / mDet;
418  mInv[0][1] = -mCD[0][1] / mDet;
419  mInv[1][1] = +mCD[0][0] / mDet;
420  sTuple.type = 2;
421  sTuple.A = (vCD[0]*mInv[0][0] + vCD[1]*mInv[1][0]) / (iPhi*iEta);
422  sTuple.B = (vCD[0]*mInv[0][1] + vCD[1]*mInv[1][1]) / (iPhi*iEta);
423  sTuple.nBins = mCD[0][0];
424  sTuple.f3 = mCD[0][1];
425  sTuple.f3sq = mCD[1][1];
426  sTuple.sig2 = vCD[0];
427  sTuple.sig2f3 = vCD[1];
428  scaleTuple->Fill(&sTuple.type);
429  BCD += sTuple.B;
430  nBCD += 1;
431  }
432  mDet = mPP[0][0]*mPP[1][1] - mPP[0][1]*mPP[1][0];
433  if (mDet > 0.00000001) {
434  mInv[0][0] = +mPP[1][1] / mDet;
435  mInv[1][0] = -mPP[1][0] / mDet;
436  mInv[0][1] = -mPP[0][1] / mDet;
437  mInv[1][1] = +mPP[0][0] / mDet;
438  sTuple.type = 3;
439  sTuple.A = (vPP[0]*mInv[0][0] + vPP[1]*mInv[1][0]) / (iPhi*iEta);
440  sTuple.B = (vPP[0]*mInv[0][1] + vPP[1]*mInv[1][1]) / (iPhi*iEta);
441  sTuple.nBins = mPP[0][0];
442  sTuple.f3 = mPP[0][1];
443  sTuple.f3sq = mPP[1][1];
444  sTuple.sig2 = vPP[0];
445  sTuple.sig2f3 = vPP[1];
446  scaleTuple->Fill(&sTuple.type);
447  BPP += sTuple.B;
448  nBPP += 1;
449  }
450  mDet = mMM[0][0]*mMM[1][1] - mMM[0][1]*mMM[1][0];
451  if (mDet > 0.00000001) {
452  mInv[0][0] = +mMM[1][1] / mDet;
453  mInv[1][0] = -mMM[1][0] / mDet;
454  mInv[0][1] = -mMM[0][1] / mDet;
455  mInv[1][1] = +mMM[0][0] / mDet;
456  sTuple.type = 4;
457  sTuple.A = (vMM[0]*mInv[0][0] + vMM[1]*mInv[1][0]) / (iPhi*iEta);
458  sTuple.B = (vMM[0]*mInv[0][1] + vMM[1]*mInv[1][1]) / (iPhi*iEta);
459  sTuple.nBins = mMM[0][0];
460  sTuple.f3 = mMM[0][1];
461  sTuple.f3sq = mMM[1][1];
462  sTuple.sig2 = vMM[0];
463  sTuple.sig2f3 = vMM[1];
464  scaleTuple->Fill(&sTuple.type);
465  BMM += sTuple.B;
466  nBMM += 1;
467  }
468  mDet = mPM[0][0]*mPM[1][1] - mPM[0][1]*mPM[1][0];
469  if (mDet > 0.00000001) {
470  mInv[0][0] = +mPM[1][1] / mDet;
471  mInv[1][0] = -mPM[1][0] / mDet;
472  mInv[0][1] = -mPM[0][1] / mDet;
473  mInv[1][1] = +mPM[0][0] / mDet;
474  sTuple.type = 5;
475  sTuple.A = (vPM[0]*mInv[0][0] + vPM[1]*mInv[1][0]) / (iPhi*iEta);
476  sTuple.B = (vPM[0]*mInv[0][1] + vPM[1]*mInv[1][1]) / (iPhi*iEta);
477  sTuple.nBins = mPM[0][0];
478  sTuple.f3 = mPM[0][1];
479  sTuple.f3sq = mPM[1][1];
480  sTuple.sig2 = vPM[0];
481  sTuple.sig2f3 = vPM[1];
482  scaleTuple->Fill(&sTuple.type);
483  BPM += sTuple.B;
484  nBPM += 1;
485  }
486  if (nSs > 0) {
487 // NSig->SetBinContent(iPhi,iEta,Ss/nSs);
488  NSigErrors->SetBinContent(iPhi,iEta,sqrt(DSs)/nDSs);
489  }
490  if (nSd > 0) {
491 // NDel->SetBinContent(iPhi,iEta,Sd/nSs);
492  NDelErrors->SetBinContent(iPhi,iEta,sqrt(DSd)/nDSd);
493  }
494  if (nSp > 0) {
495 // NPlus->SetBinContent(iPhi,iEta,Sp/nSp);
496  NPlusErrors->SetBinContent(iPhi,iEta,sqrt(DSp)/nDSp);
497  }
498  if (nSm > 0) {
499 // NMinus->SetBinContent(iPhi,iEta,Sm/nSm);
500  NMinusErrors->SetBinContent(iPhi,iEta,sqrt(DSm)/nDSm);
501  }
502  if (nSc > 0) {
503 // NPlusMinus->SetBinContent(iPhi,iEta,Sc/nSc);
504  NPlusMinusErrors->SetBinContent(iPhi,iEta,sqrt(DSc)/nDSc);
505  }
506  }
507  }
508  STuple.type = 1;
509  STuple.B = BCI;
510  STuple.nBins = nBCI;
511  sumTuple->Fill(&STuple.type);
512  STuple.type = 2;
513  STuple.B = BCD;
514  STuple.nBins = nBCD;
515  sumTuple->Fill(&STuple.type);
516  STuple.type = 3;
517  STuple.B = BPP;
518  STuple.nBins = nBPP;
519  sumTuple->Fill(&STuple.type);
520  STuple.type = 4;
521  STuple.B = BMM;
522  STuple.nBins = nBMM;
523  sumTuple->Fill(&STuple.type);
524  STuple.type = 5;
525  STuple.B = BPM;
526  STuple.nBins = nBPM;
527  sumTuple->Fill(&STuple.type);
528 
529  BCI = BCI / nBCI;
530  BCD = BCD / nBCD;
531  BPP = BPP / nBPP;
532  BMM = BMM / nBMM;
533  BPM = BPM / nBPM;
534  double delS, corr, a, b, f;
535  double fFull = (pow(mEtaMax,3)-pow(mEtaMin,3))/(3*(mEtaMax-mEtaMin));
536  for (int jPhi=0;jPhi<mNPhiBins;jPhi++) {
537  int iPhi = jPhi + 1;
538  for (int jEta=0;jEta<mNEtaBins;jEta++) {
539  int iEta = jEta + 1;
540  if (vecCI[jPhi][jEta][0] > 0) {
541  a = vecCI[jPhi][jEta][1]/vecCI[jPhi][jEta][0];
542  b = BCI * iPhi*iEta;
543  f = vecCI[jPhi][jEta][2]/vecCI[jPhi][jEta][0];
544 // delS = a + b*(1.0/3 - f);
545  corr = b*(f - fFull);
546 // corr = b*f;
547  delS = a - corr;
548  NSig->SetBinContent(iPhi,iEta,delS);
549  NSigCorrection->SetBinContent(iPhi,iEta,corr);
550  } else {
551  NSig->SetBinContent(iPhi,iEta,0);
552  }
553  if (vecCD[jPhi][jEta][0] > 0) {
554  a = vecCD[jPhi][jEta][1]/vecCD[jPhi][jEta][0];
555  b = BCD * iPhi*iEta;
556  f = vecCD[jPhi][jEta][2]/vecCD[jPhi][jEta][0];
557 // delS = a + b*(1.0/3 - f);
558  corr = b*(f - fFull);
559 // corr = b*f;
560  delS = a - corr;
561  NDel->SetBinContent(iPhi,iEta,delS);
562  NDelCorrection->SetBinContent(iPhi,iEta,corr);
563  } else {
564  NDel->SetBinContent(iPhi,iEta,0);
565  }
566  if (vecPP[jPhi][jEta][0] > 0) {
567  a = vecPP[jPhi][jEta][1]/vecPP[jPhi][jEta][0];
568  b = BPP * iPhi*iEta;
569  f = vecPP[jPhi][jEta][2]/vecPP[jPhi][jEta][0];
570 // delS = a + b*(1.0/3 - f);
571  corr = b*(f - fFull);
572 // corr = b*f;
573  delS = a - corr;
574  NPlus->SetBinContent(iPhi,iEta,delS);
575  NPlusCorrection->SetBinContent(iPhi,iEta,corr);
576  } else {
577  NPlus->SetBinContent(iPhi,iEta,0);
578  }
579  if (vecMM[jPhi][jEta][0] > 0) {
580  a = vecMM[jPhi][jEta][1]/vecMM[jPhi][jEta][0];
581  b = BMM * iPhi*iEta;
582  f = vecMM[jPhi][jEta][2]/vecMM[jPhi][jEta][0];
583 // delS = a + b*(1.0/3 - f);
584  corr = b*(f - fFull);
585 // corr = b*f;
586  delS = a - corr;
587  NMinus->SetBinContent(iPhi,iEta,delS);
588  NMinusCorrection->SetBinContent(iPhi,iEta,corr);
589  } else {
590  NMinus->SetBinContent(iPhi,iEta,0);
591  }
592  if (vecPM[jPhi][jEta][0] > 0) {
593  a = vecPM[jPhi][jEta][1]/vecPM[jPhi][jEta][0];
594  b = BPM * iPhi*iEta;
595  f = vecPM[jPhi][jEta][2]/vecPM[jPhi][jEta][0];
596 // delS = a + b*(1.0/3 - f);
597  corr = b*(f - fFull);
598 // corr = b*f;
599  delS = a - corr;
600  NPlusMinus->SetBinContent(iPhi,iEta,delS);
601  NPlusMinusCorrection->SetBinContent(iPhi,iEta,corr);
602  } else {
603  NPlusMinus->SetBinContent(iPhi,iEta,0);
604  }
605  }
606  }
607 }
608 void StEStructSigmas::PHistograms() {
609  char buffer[255];
610  binTupleStruct bTuple;
611 
612  TH2F *hnBins = (TH2F *) gDirectory->Get("nBins");
613  TH2F *hoffset = (TH2F *) gDirectory->Get("offset");
614  TH2F *hfUnique = (TH2F *) gDirectory->Get("fUnique");
615 
616  // We have accumulated over events.
617  // For every bin at each scale we have summed up all the small bins.
618  // Now accumulate sum, n, n^2 for all bins.
619 
620  for (int iType=0;iType<3;iType++) {
621  PSig[iType]->Reset();
622  PDel[iType]->Reset();
623  PPlus[iType]->Reset();
624  PMinus[iType]->Reset();
625  PPlusMinus[iType]->Reset();
626  }
627  PSig[3]->Reset();
628 
629  PSigErrors->Reset();
630  PDelErrors->Reset();
631  PPlusErrors->Reset();
632  PMinusErrors->Reset();
633  PPlusMinusErrors->Reset();
634 
635  SPtHat->Reset();
636  PPtHat->Reset();
637  MPtHat->Reset();
638  sigSPtHat->Reset();
639  sigPPtHat->Reset();
640  sigMPtHat->Reset();
641 
642  sigSPtHatErrors->Reset();
643  sigPPtHatErrors->Reset();
644  sigMPtHatErrors->Reset();
645 
646 
647  TH1D *hTotEvents[6];
648 
649  TH1D *hNSum[16];
650  TH1D *hNDiff[2];
651  TH1D *hNPlus[5];
652  TH1D *hNMinus[5];
653  TH1D *hNPlusMinus[8];
654 
655  TH1D *hPSum[21];
656  TH1D *hPDiff[16];
657  TH1D *hPPlus[11];
658  TH1D *hPMinus[11];
659  TH1D *hPPlusMinus[17];
660 
661  sprintf(buffer,"%sTotalEvents_%s", mpreFix, mKey);
662  hTotEvents[0] = (TH1D *) gDirectory->Get(buffer);
663  sprintf(buffer,"%sTotalSumEvents_%s", mpreFix, mKey);
664  hTotEvents[1] = (TH1D *) gDirectory->Get(buffer);
665  sprintf(buffer,"%sTotalPlusEvents_%s", mpreFix, mKey);
666  hTotEvents[2] = (TH1D *) gDirectory->Get(buffer);
667  sprintf(buffer,"%sTotalMinusEvents_%s", mpreFix, mKey);
668  hTotEvents[3] = (TH1D *) gDirectory->Get(buffer);
669  sprintf(buffer,"%sTotalPlusMinusEvents_%s", mpreFix, mKey);
670  hTotEvents[4] = (TH1D *) gDirectory->Get(buffer);
671  sprintf(buffer, "%sTotalEventsForSig2Dyn_%s", mpreFix, mKey);
672  hTotEvents[5] = (TH1D *) gDirectory->Get(buffer);
673 
674 
675  for (int jStat=0;jStat<16;jStat++) {
676  sprintf( buffer,"%sNSum%s_%i", mpreFix, mKey, jStat );
677  hNSum[jStat] = (TH1D *) gDirectory->Get(buffer);
678  }
679  for (int jStat=0;jStat<2;jStat++) {
680  sprintf( buffer,"%sNDiff%s_%i", mpreFix, mKey, jStat );
681  hNDiff[jStat] = (TH1D *) gDirectory->Get(buffer);
682  }
683  for (int jStat=0;jStat<5;jStat++) {
684  sprintf( buffer,"%sNPlus%s_%i", mpreFix, mKey, jStat );
685  hNPlus[jStat] = (TH1D *) gDirectory->Get(buffer);
686  }
687  for (int jStat=0;jStat<5;jStat++) {
688  sprintf( buffer,"%sNMinus%s_%i", mpreFix, mKey, jStat );
689  hNMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
690  }
691  for (int jStat=0;jStat<8;jStat++) {
692  sprintf( buffer,"%sNPlusMinus%s_%i", mpreFix, mKey, jStat );
693  hNPlusMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
694  }
695 
696 
697  for (int jStat=0;jStat<21;jStat++) {
698  sprintf(buffer,"%sPSum%s_%i", mpreFix, mKey, jStat);
699  hPSum[jStat] = (TH1D *) gDirectory->Get(buffer);
700  }
701  for (int jStat=0;jStat<16;jStat++) {
702  sprintf(buffer,"%sPDiff%s_%i", mpreFix, mKey, jStat);
703  hPDiff[jStat] = (TH1D *) gDirectory->Get(buffer);
704  }
705  for (int jStat=0;jStat<11;jStat++) {
706  sprintf(buffer,"%sPPlus%s_%i", mpreFix, mKey, jStat);
707  hPPlus[jStat] = (TH1D *) gDirectory->Get(buffer);
708  }
709  for (int jStat=0;jStat<11;jStat++) {
710  sprintf(buffer,"%sPMinus%s_%i", mpreFix, mKey, jStat);
711  hPMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
712  }
713  for (int jStat=0;jStat<17;jStat++) {
714  sprintf(buffer,"%sPPlusMinus%s_%i", mpreFix, mKey, jStat);
715  hPPlusMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
716  }
717 
718  int mBin, oBin;
719  int nBad = 0;
720 
721  // Currently working with `small' event samples and calculating \hat{p_t}
722  // for small bins of peripheral events is problematic. Since we are dealing
723  // with Hijing I try using an inclusive \hat{p_t}
724  // Maybe bad for data where mean p_t may depend on detector location.
725 
726 
727  int iBinFS = (int) (hnBins->GetBinContent(mNPhiBins,mNEtaBins)+hoffset->GetBinContent(mNPhiBins,mNEtaBins));
728  double nplus = (double) hNSum[0]->GetBinContent(iBinFS);
729  double nminus = (double) hNSum[1]->GetBinContent(iBinFS);
730  double pplus = (double) hPSum[1]->GetBinContent(iBinFS);
731  double pminus = (double) hPSum[2]->GetBinContent(iBinFS);
732  double psq = (double) hPSum[0]->GetBinContent(iBinFS);
733  double psig = (double) hPSum[17]->GetBinContent(iBinFS);
734  double nsig = (double) hPSum[16]->GetBinContent(iBinFS);
735  double psigsq = (double) hPSum[18]->GetBinContent(iBinFS);
736  double nsum = nplus + nminus;
737  double psum = pplus + pminus;
738 
739  for (int iPhi=1;iPhi<=mNPhiBins;iPhi++) {
740  for (int iEta=1;iEta<=mNEtaBins;iEta++) {
741  double NS[16], ND[2], NP[5], NM[5], NC[8];
742  double PS[21], PD[16], PP[11], PM[11], PC[17];
743 
744  double pHat[3], pHat2[3], pHatSig2, sigHatSig2;
745  double sigSq0, sigSq1, sigSq2, sigHat, DsigHat;
746  double Ss[] = {0, 0, 0, 0}, nSs = 0, nS1 = 0, nS3 = 0, DSs = 0, nDSs = 0, hSs = 0, hDSs = 0;
747  double Sd[] = {0, 0, 0}, nSd = 0, DSd = 0, nDSd = 0;
748  double Sp[] = {0, 0, 0}, nSp = 0, DSp = 0, nDSp = 0, hSp = 0, hDSp = 0;
749  double Sm[] = {0, 0, 0}, nSm = 0, DSm = 0, nDSm = 0, hSm = 0, hDSm = 0;
750  double Sc[] = {0, 0, 0}, nSc = 0, nDSc = 0;
751 
752  double hatSs = 0;
753  double hatSp = 0;
754  double hatSm = 0;
755 
756  double totEvents, sumEvents, sig2Events, plusEvents, minusEvents, plusMinusEvents;
757 
758  mBin = (int) hnBins->GetBinContent(iPhi,iEta);
759  oBin = (int) hoffset->GetBinContent(iPhi,iEta);
760 
761  int iEtaBin = 0, iPhiBin = 1;
762  int nEtaBin = getNumEtaBins(iEta);
763  float phi1, phi2, eta1, eta2;
764  double nTracksTotal = 0, nBinsTotal = 0;
765  for (int iBin=oBin+1;iBin<=oBin+mBin;iBin++) {
766  iEtaBin++;
767  if (iEtaBin > nEtaBin) {
768  iEtaBin = 1;
769  iPhiBin++;
770  }
771  phi1 = 0 + 2*3.1415926*getPhiStart(iPhiBin,iPhi)/mNPhiBins;
772  phi2 = phi1 + iPhi * 2*3.1415926 / mNPhiBins;
773  eta1 = mEtaMin + (mEtaMax-mEtaMin)*getEtaStart(iEtaBin,iEta)/mNEtaBins;
774  eta2 = eta1 + iEta * (mEtaMax-mEtaMin) / mNEtaBins;
775  totEvents = hTotEvents[0]->GetBinContent(iBin);
776  sumEvents = hTotEvents[1]->GetBinContent(iBin);
777  sig2Events = hTotEvents[5]->GetBinContent(iBin);
778  if (sumEvents > 0) {
779  for (int jStat=0;jStat<16;jStat++) {
780  NS[jStat] = hNSum[jStat]->GetBinContent(iBin);
781  }
782  for (int jStat=0;jStat<2;jStat++) {
783  ND[jStat] = hNDiff[jStat]->GetBinContent(iBin);
784  }
785  for (int jStat=0;jStat<21;jStat++) {
786  PS[jStat] = hPSum[jStat]->GetBinContent(iBin);
787  }
788  for (int jStat=0;jStat<16;jStat++) {
789  PD[jStat] = hPDiff[jStat]->GetBinContent(iBin);
790  }
791  double NSumSig2 = PS[16];
792  double NPlus = NS[0];
793  double PPlus = PS[1];
794  double NMinus = NS[1];
795  double PMinus = PS[2];
796  double NSum = NPlus + NMinus;
797  double PSum = PPlus + PMinus;
798  nTracksTotal += NSum;
799  nBinsTotal += totEvents;
800  if (NSum > 0) {
801  // Lower case values are for largest bin at this centrality.
802  // Upper case are for the current bin.
803  if (nsum > 0) {
804 // pHat[0] = PSum / NSum;
805  pHat[0] = psum / nsum;
806  pHat2[0] = pHat[0]*pHat[0];
807  if (nplus > 0) {
808 // pHat[1] = PPlus / NPlus;
809  pHat[1] = pplus / nplus;
810  } else {
811  pHat[1] = 0;
812  }
813  pHat2[1] = pHat[1]*pHat[1];
814  if (nminus > 0) {
815 // pHat[2] = PMinus / NMinus;
816  pHat[2] = pminus / nminus;
817  } else {
818  pHat[2] = 0;
819  }
820  pHat2[2] = pHat[2]*pHat[2];
821  // Using all tracks for this centrality in calculating \sigma^2_{\hat{p_t}}
822  // is fine when we sum over all bins.
823  sigHat = psq/nsum - (psum/nsum)*(psum/nsum);
824 
825 // if (NSig > 0) {
826 // pHatSig2 = PSig / NSig;
827 // } else {
828 // pHatSig2 = 0;
829 // }
830  if (nsig > 0) {
831  pHatSig2 = psig / nsig;
832  sigHatSig2 = psigsq / nsig - (psig/nsig)*(psig/nsig);
833  } else {
834  pHatSig2 = 0;
835  sigHatSig2 = 0;
836  }
837  }
838 
839 
840  hatSs += pHat[0];
841  hSs += sigHat;
842  DsigHat = 4*pHat2[0]*sigHat / (NSum*sumEvents);
843  hDSs += DsigHat;
844  sigSq0 = PS[8] - 2*pHat[1]*PS[6] - 2*pHat[2]*PS[7]
845  + pHat2[1]*NS[5] + 2*pHat[1]*pHat[2]*NS[6] + pHat2[2]*NS[7];
846  Ss[0] += sigSq0 - sigHat;
847  sigSq1 = PS[5] - 2*pHat[1]*PS[3] - 2*pHat[2]*PS[4]
848  + pHat2[1]*NS[2] + 2*pHat[1]*pHat[2]*NS[3] + pHat2[2]*NS[4];
849  // Create \Delta^\sigma^2 / \bar n to compare with <dpt dpt>
850  if (NSum > 1) {
851  Ss[1] += sigSq1/NSum - sigHat;
852 // Ss[1] += (sigSq1/NSum - sigHat) / (NSum/totEvents);
853  nS1 += 1;
854  } else {
855  nBad += 1;
856  }
857 // Doing some tests for Jamie to alleive his concerns over nu_dynamical calculation.
858 // My original comparison in 2004 was reasonable, but in comparing to to EbyE
859 // energy dependence paper I should have used their newer definition, which
860 // here I am calling 2005.
861  if (NSumSig2 > 0) { // Here is calculation with old reference.
862  sigSq2 = PS[11] - 2*pHatSig2*PS[10] + pHatSig2*pHatSig2*PS[9];
863  Ss[2] += (NSumSig2 -1)* (sigSq2 - sigHatSig2*NS[15]);
864  }
865  if (NSumSig2 > 1) { // Here is calculation with 2005 reference.
866  // Create <dpt dpt> histograms with or without \bar{n} normalization.
867  // Choose global or local pHat when calculating pHatSig2 above.
868  sigSq2 = PS[11] - 2*pHatSig2*PS[10] + pHatSig2*pHatSig2*PS[9];
869 // Ss[3] += (NSum/totEvents)*(sigSq2 - (PS[19] -2*pHatSig2*PS[20] + pHatSig2*pHatSig2*NS[15]))/sig2Events;
870  Ss[3] += (sigSq2 - (PS[19] -2*pHatSig2*PS[20] + pHatSig2*pHatSig2*NS[15]))/sig2Events;
871  nS3 += 1;
872  }
873  nSs += 1;
874  bTuple.type = 6;
875  bTuple.phiScale = iPhi;
876  bTuple.etaScale = iEta;
877  bTuple.phi = (phi1+phi2)/2;
878  bTuple.eta = (eta1+eta2)/2;
879  bTuple.sig2 = sigSq0 - sigHat;
880  bTuple.sig2_1 = sigSq1/NSum - sigHat;
881  bTuple.sig2_2 = NSum * (sigSq2 - sigHat*NS[15]);
882  bTuple.nbar = NSum;
883  bTuple.events = sumEvents;
884  bTuple.f3 = 0;
885  binTuple->Fill(&bTuple.type);
886 
887  DSs += (PS[13] - 4*pHat[0]*PS[12] + 6*pHat2[0]*PS[8]
888  - 4*pHat[0]*pHat2[0]*PSum + pHat2[0]*pHat2[0]*NSum) *NSum/sumEvents;
889  nDSs += sqrt(hfUnique->GetBinContent(iPhi,iEta));
890 
891  sigSq0 = PD[5] - 2*pHat[1]*PD[3] + 2*pHat[2]*PD[4]
892  + pHat2[1]*NS[5] - 2*pHat[1]*pHat[2]*NS[6] + pHat2[2]*NS[7];
893  Sd[0] += sigSq0 - sigHat;
894  sigSq1 = PD[2] - 2*pHat[1]*PD[0] + 2*pHat[2]*PD[1]
895  + pHat2[1]*NS[2] - 2*pHat[1]*pHat[2]*NS[3] + pHat2[2]*NS[4];
896  Sd[1] += sigSq1/NSum - sigHat;
897  sigSq2 = PD[8] - 2*pHat[1]*PD[6] + 2*pHat[2]*PD[7]
898  + pHat2[1]*NS[8] - 2*pHat[1]*pHat[2]*NS[9] + pHat2[2]*NS[10];
899  Sd[2] += NSum * (sigSq2 - sigHat*NS[15]);
900  nSd += 1;
901 // Ss[0] += PD[2]/NSum;
902 // Ss[1] += 2*pHat[1]*PD[0]/NSum;
903 // Ss[2] += 2*pHat[2]*PD[1]/NSum;
904 // Sd[0] += pHat2[1]*NS[2]/NSum;
905 // Sd[1] += 2*pHat[1]*pHat[2]*NS[3]/NSum;
906 // Sd[2] += pHat2[2]*NS[4]/NSum;
907 
908  DSd += (PD[11] - 4*pHat[0]*PD[12] + 6*pHat2[0]*PD[13]
909  - 4*pHat[0]*pHat2[0]*PD[14] + pHat2[0]*pHat2[0]*PD[15]) *NSum/sumEvents;
910  nDSd += sqrt(hfUnique->GetBinContent(iPhi,iEta));
911  bTuple.type = 7;
912  bTuple.phiScale = iPhi;
913  bTuple.etaScale = iEta;
914  bTuple.phi = (phi1+phi2)/2;
915  bTuple.eta = (eta1+eta2)/2;
916  bTuple.sig2 = sigSq0 - sigHat;
917  bTuple.sig2_1 = sigSq1/NSum - sigHat;
918  bTuple.sig2_2 = NSum * (sigSq2 - sigHat*NS[15]);
919  bTuple.nbar = NSum;
920  bTuple.events = sumEvents;
921  bTuple.f3 = 0;
922  binTuple->Fill(&bTuple.type);
923  }
924  }
925 
926  plusEvents = hTotEvents[2]->GetBinContent(iBin);
927  pHat[1] = 0;
928  if (plusEvents > 0) {
929  for (int jStat=0;jStat<5;jStat++) {
930  NP[jStat] = hNPlus[jStat]->GetBinContent(iBin) / plusEvents;
931  }
932  for (int jStat=0;jStat<11;jStat++) {
933  PP[jStat] = hPPlus[jStat]->GetBinContent(iBin) / plusEvents;
934  }
935  if (NP[0] > 0) {
936  pHat[1] = PP[1] / NP[0];
937  pHat2[1] = pHat[1]*pHat[1];
938  hatSp += pHat[1];
939  sigHat = PP[0]/NP[0] - pHat[1]*pHat[1];
940  hSp += sigHat;
941  DsigHat = 4*pHat2[1]*sigHat / (NP[0]*plusEvents);
942  hDSp += DsigHat;
943  sigSq0 = PP[5] - PP[1]*PP[1]/NP[0];
944  Sp[0] += sigSq0 - sigHat;
945  sigSq1 = PP[4] - 2*pHat[1]*PP[2] + pHat2[1]*NP[1];
946  Sp[1] += sigSq1/NP[0] - sigHat;
947  sigSq2 = PP[6] - 2*pHat[1]*PP[3] + pHat2[1];
948  Sp[2] += NP[0] * (sigSq2 - sigHat*NP[4]);
949  nSp += 1;
950 
951  DSp += (PP[8] - 4*pHat[1]*PP[7] + 6*pHat2[1]*PP[5]
952  - 4*pHat[1]*pHat2[1]*PP[1] + pHat2[1]*pHat2[1]*NP[0]) *NP[0]/plusEvents;
953  nDSp += sqrt(hfUnique->GetBinContent(iPhi,iEta));
954  bTuple.type = 8;
955  bTuple.phiScale = iPhi;
956  bTuple.etaScale = iEta;
957  bTuple.phi = (phi1+phi2)/2;
958  bTuple.eta = (eta1+eta2)/2;
959  bTuple.sig2 = sigSq0 - sigHat;
960  bTuple.sig2_1 = sigSq1/NP[0] - sigHat;
961  bTuple.sig2_2 = NP[0] * (sigSq2 - sigHat*NP[4]);
962  bTuple.nbar = NP[0];
963  bTuple.events = plusEvents;
964  bTuple.f3 = 0;
965  binTuple->Fill(&bTuple.type);
966  }
967  } else {
968  for (int jStat=0;jStat<5;jStat++) {
969  NP[jStat] = 0;
970  }
971  for (int jStat=0;jStat<11;jStat++) {
972  PP[jStat] = 0;
973  }
974  }
975 
976  minusEvents = hTotEvents[3]->GetBinContent(iBin);
977  pHat[2] = 0;
978  if (minusEvents > 0) {
979  for (int jStat=0;jStat<5;jStat++) {
980  NM[jStat] = hNMinus[jStat]->GetBinContent(iBin) / minusEvents;
981  }
982  for (int jStat=0;jStat<11;jStat++) {
983  PM[jStat] = hPMinus[jStat]->GetBinContent(iBin) / minusEvents;
984  }
985  if (NM[0] > 0) {
986  pHat[2] = PM[1] / NM[0];
987  pHat2[2] = pHat[2]*pHat[2];
988  hatSm += pHat[2];
989  sigHat = PM[0]/NM[0] - pHat[2]*pHat[2];
990  hSm += sigHat;
991  DsigHat = 4*pHat2[2]*sigHat / (NM[0]*minusEvents);
992  hDSm += DsigHat;
993  sigSq0 = PM[5] - PM[1]*PM[1]/NM[0];
994  Sm[0] += sigSq0 - sigHat;
995  sigSq1 = PM[4] - 2*pHat[2]*PM[2] + pHat2[2]*NM[1];
996  Sm[1] += sigSq1/NM[0] - sigHat;
997  sigSq2 = PM[6] - 2*pHat[2]*PM[3] + pHat2[2];
998  Sm[2] += NM[0] * (sigSq2 - sigHat*NM[4]);
999  nSm += 1;
1000 
1001  DSm += (PM[8] - 4*pHat[2]*PM[7] + 6*pHat2[2]*PM[5]
1002  - 4*pHat[2]*pHat2[2]*PM[1] + pHat2[2]*pHat2[2]*NM[0]) * NM[0] / minusEvents;
1003  nDSm += sqrt(hfUnique->GetBinContent(iPhi,iEta));
1004  bTuple.type = 9;
1005  bTuple.phiScale = iPhi;
1006  bTuple.etaScale = iEta;
1007  bTuple.phi = (phi1+phi2)/2;
1008  bTuple.eta = (eta1+eta2)/2;
1009  bTuple.sig2 = sigSq0 - sigHat;
1010  bTuple.sig2_1 = sigSq1/NM[0] - sigHat;
1011  bTuple.sig2_2 = NM[0] * (sigSq2 - sigHat*NM[4]);
1012  bTuple.nbar = NM[0];
1013  bTuple.events = minusEvents;
1014  bTuple.f3 = 0;
1015  binTuple->Fill(&bTuple.type);
1016  }
1017  } else {
1018  for (int jStat=0;jStat<5;jStat++) {
1019  NM[jStat] = 0;
1020  }
1021  for (int jStat=0;jStat<11;jStat++) {
1022  PM[jStat] = 0;
1023  }
1024  }
1025 
1026 // Accumulation of hNPlusMinus required (n+ > 0) && (n- > 0)
1027 // I think this is wrong. Instead if I require (n+ > 0) || (n- > 0)
1028 // I can retrieve appropriate terms from NP and NM, at least for the
1029 // phiPt case.
1030  plusMinusEvents = hTotEvents[1]->GetBinContent(iBin);
1031  if (plusMinusEvents > 0) {
1032  for (int jStat=0;jStat<8;jStat++) {
1033  NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1034  }
1035  for (int jStat=0;jStat<17;jStat++) {
1036  PC[jStat] = hPPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1037  }
1038 
1039  sigSq0 = PC[16] - pHat[2]*PC[10] - pHat[1]*PC[11] + pHat[1]*pHat[2]*NC[5];
1040  Sc[0] += sigSq0;
1041  sigSq1 = 0;
1042  if (NP[0]*NM[0] > 0) {
1043  sigSq1 = PC[14] - pHat[2]*PC[5] - pHat[1]*PC[4] + pHat[1]*pHat[2]*NC[2];
1044  Sc[1] += sigSq1/sqrt(NP[0]*NM[0]);
1045  }
1046  sigSq2 = PC[15] - pHat[2]*PC[2] - pHat[1]*PC[3] + pHat[1]*pHat[2];
1047  Sc[2] += sigSq2*sqrt(NP[0]*NM[0]);
1048  nSc += 1;
1049  nDSc += sqrt(hfUnique->GetBinContent(iPhi,iEta));
1050  bTuple.type = 10;
1051  bTuple.phiScale = iPhi;
1052  bTuple.etaScale = iEta;
1053  bTuple.phi = (phi1+phi2)/2;
1054  bTuple.eta = (eta1+eta2)/2;
1055  bTuple.sig2 = sigSq0;
1056  bTuple.sig2_1 = sigSq1/sqrt(NP[0]*NM[0]);
1057  bTuple.sig2_2 = sigSq2*sqrt(NP[0]*NM[0]);
1058  bTuple.nbar = sqrt(NP[0]*NM[0]);
1059  bTuple.events = plusMinusEvents;
1060  bTuple.f3 = 0;
1061  binTuple->Fill(&bTuple.type);
1062  }
1063 // plusMinusEvents = hTotEvents[1]->GetBinContent(iBin);
1064 // if (plusMinusEvents > 0) {
1065 // for (int jStat=0;jStat<8;jStat++) {
1066 // NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1067 // }
1068 // for (int jStat=0;jStat<17;jStat++) {
1069 // PC[jStat] = hPPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1070 // }
1071 // if (NC[0] > 0) {
1072 // pHat[1] = PC[0] / NC[0];
1073 // } else {
1074 // pHat[1] = 0;
1075 // }
1076 // if (NC[1] > 0) {
1077 // pHat[2] = PC[1] / NC[1];
1078 // } else {
1079 // pHat[2] = 0;
1080 // }
1081 // sigSq = PC[16] - pHat[2]*PC[10] - pHat[1]*PC[11] + pHat[1]*pHat[2]*NC[5];
1082 // Sc[0] += sigSq;
1083 // sigSq = PC[14] - pHat[2]*PC[5] - pHat[1]*PC[4] + pHat[1]*pHat[2]*NC[2];
1084 // Sc[1] += sigSq/sqrt(NC[0]*NC[1]);
1085 // sigSq = PC[15] - pHat[2]*PC[2] - pHat[1]*PC[3] + pHat[1]*pHat[2];
1086 // Sc[2] += sigSq*sqrt(NC[0]*NC[1]);
1087 // nSc += 1;
1088 // nDSc += sqrt(hfUnique->GetBinContent(iPhi,iEta));
1089 // }
1090  }
1091 
1092  if (nSs > 0) {
1093  PSig[0]->SetBinContent(iPhi,iEta,Ss[0]/nSs);
1094  PSigErrors->SetBinContent(iPhi,iEta,sqrt(DSs)/nDSs);
1095  SPtHat->SetBinContent(iPhi,iEta,hatSs/nSs);
1096  sigSPtHat->SetBinContent(iPhi,iEta,hSs/nSs);
1097  sigSPtHatErrors->SetBinContent(iPhi,iEta,sqrt(hDSs)/nDSs);
1098  }
1099  if (nS1 > 0) {
1100  PSig[1]->SetBinContent(iPhi,iEta,Ss[1]/nS1);
1101 // PSig[1]->SetBinContent(iPhi,iEta,(Ss[1]/nS1)*1.0/(nTracksTotal/nBinsTotal));
1102  }
1103  if (nS3 > 0) {
1104  PSig[2]->SetBinContent(iPhi,iEta,Ss[2]/nSs);
1105  PSig[3]->SetBinContent(iPhi,iEta,Ss[3]/nS3);
1106  }
1107  if (nSd > 0) {
1108  PDel[0]->SetBinContent(iPhi,iEta,Sd[0]/nSd);
1109  PDel[1]->SetBinContent(iPhi,iEta,Sd[1]/nSd);
1110  PDel[2]->SetBinContent(iPhi,iEta,Sd[2]/nSd);
1111  PDelErrors->SetBinContent(iPhi,iEta,sqrt(DSd)/nDSd);
1112  }
1113 
1114  if (nSp > 0) {
1115  PPlus[0]->SetBinContent(iPhi,iEta,Sp[0]/nSp);
1116  PPlus[1]->SetBinContent(iPhi,iEta,Sp[1]/nSp);
1117  PPlus[2]->SetBinContent(iPhi,iEta,Sp[2]/nSp);
1118  PPlusErrors->SetBinContent(iPhi,iEta,sqrt(DSp)/nDSp);
1119  PPtHat->SetBinContent(iPhi,iEta,hatSp/nSp);
1120  sigPPtHat->SetBinContent(iPhi,iEta,hSp/nSp);
1121  sigPPtHatErrors->SetBinContent(iPhi,iEta,sqrt(hDSp)/nDSp);
1122  }
1123 
1124  if (nSm > 0) {
1125  PMinus[0]->SetBinContent(iPhi,iEta,Sm[0]/nSm);
1126  PMinus[1]->SetBinContent(iPhi,iEta,Sm[1]/nSm);
1127  PMinus[2]->SetBinContent(iPhi,iEta,Sm[2]/nSm);
1128  PMinusErrors->SetBinContent(iPhi,iEta,sqrt(DSm)/nDSm);
1129  MPtHat->SetBinContent(iPhi,iEta,hatSm/nSm);
1130  sigMPtHat->SetBinContent(iPhi,iEta,hSm/nSm);
1131  sigMPtHatErrors->SetBinContent(iPhi,iEta,sqrt(hDSm)/nDSm);
1132  }
1133 
1134  if (nSc > 0) {
1135  PPlusMinus[0]->SetBinContent(iPhi,iEta,Sc[0]/nSc);
1136  PPlusMinus[1]->SetBinContent(iPhi,iEta,Sc[1]/nSc);
1137  PPlusMinus[2]->SetBinContent(iPhi,iEta,Sc[2]/nSc);
1138 // Haven't found good approximation for error on +- covariance.
1139 // Use average of errors on + and - to have something.
1140  PPlusMinusErrors->SetBinContent(iPhi,iEta,sqrt((DSp+DSm)/2)/nDSs);
1141  }
1142  }
1143  }
1144 }
1145 void StEStructSigmas::PNHistograms() {
1146  char buffer[255];
1147  binTupleStruct bTuple;
1148 
1149  TH2F *hnBins = (TH2F *) gDirectory->Get("nBins");
1150  TH2F *hoffset = (TH2F *) gDirectory->Get("offset");
1151 
1152  // We have accumulated over events.
1153  // For every bin at each scale we have summed up all the small bins.
1154  // Now accumulate sum, n, n^2 for all bins.
1155 
1156  for (int iType=0;iType<3;iType++) {
1157  PNSig[iType]->Reset();
1158  PNDel[iType]->Reset();
1159  PNPlus[iType]->Reset();
1160  PNMinus[iType]->Reset();
1161  PNPlusMinus[iType]->Reset();
1162  PNMinusPlus[iType]->Reset();
1163  }
1164 
1165  PNSigErrors->Reset();
1166  PNDelErrors->Reset();
1167  PNPlusErrors->Reset();
1168  PNMinusErrors->Reset();
1169  PNPlusMinusErrors->Reset();
1170  PNMinusPlusErrors->Reset();
1171 
1172  TH1D *hTotEvents[5];
1173 
1174  TH1D *hNSum[16];
1175  TH1D *hNDiff[2];
1176  TH1D *hNPlus[5];
1177  TH1D *hNMinus[5];
1178  TH1D *hNPlusMinus[8];
1179 
1180  TH1D *hPSum[21];
1181  TH1D *hPDiff[16];
1182  TH1D *hPPlus[11];
1183  TH1D *hPMinus[11];
1184  TH1D *hPPlusMinus[17];
1185 
1186  sprintf(buffer,"%sTotalEvents_%s", mpreFix, mKey);
1187  hTotEvents[0] = (TH1D *) gDirectory->Get(buffer);
1188  sprintf(buffer,"%sTotalSumEvents_%s", mpreFix, mKey);
1189  hTotEvents[1] = (TH1D *) gDirectory->Get(buffer);
1190  sprintf(buffer,"%sTotalPlusEvents_%s", mpreFix, mKey);
1191  hTotEvents[2] = (TH1D *) gDirectory->Get(buffer);
1192  sprintf(buffer,"%sTotalMinusEvents_%s", mpreFix, mKey);
1193  hTotEvents[3] = (TH1D *) gDirectory->Get(buffer);
1194  sprintf(buffer,"%sTotalPlusMinusEvents_%s", mpreFix, mKey);
1195  hTotEvents[4] = (TH1D *) gDirectory->Get(buffer);
1196 
1197  for (int jStat=0;jStat<16;jStat++) {
1198  sprintf( buffer,"%sNSum%s_%i", mpreFix, mKey, jStat );
1199  hNSum[jStat] = (TH1D *) gDirectory->Get(buffer);
1200  }
1201  for (int jStat=0;jStat<2;jStat++) {
1202  sprintf( buffer,"%sNDiff%s_%i", mpreFix, mKey, jStat );
1203  hNDiff[jStat] = (TH1D *) gDirectory->Get(buffer);
1204  }
1205  for (int jStat=0;jStat<5;jStat++) {
1206  sprintf( buffer,"%sNPlus%s_%i", mpreFix, mKey, jStat );
1207  hNPlus[jStat] = (TH1D *) gDirectory->Get(buffer);
1208  }
1209  for (int jStat=0;jStat<5;jStat++) {
1210  sprintf( buffer,"%sNMinus%s_%i", mpreFix, mKey, jStat );
1211  hNMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
1212  }
1213  for (int jStat=0;jStat<8;jStat++) {
1214  sprintf( buffer,"%sNPlusMinus%s_%i", mpreFix, mKey, jStat );
1215  hNPlusMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
1216  }
1217 
1218 
1219  for (int jStat=0;jStat<21;jStat++) {
1220  sprintf(buffer,"%sPSum%s_%i", mpreFix, mKey, jStat);
1221  hPSum[jStat] = (TH1D *) gDirectory->Get(buffer);
1222  }
1223  for (int jStat=0;jStat<16;jStat++) {
1224  sprintf(buffer,"%sPDiff%s_%i", mpreFix, mKey, jStat);
1225  hPDiff[jStat] = (TH1D *) gDirectory->Get(buffer);
1226  }
1227  for (int jStat=0;jStat<11;jStat++) {
1228  sprintf(buffer,"%sPPlus%s_%i", mpreFix, mKey, jStat);
1229  hPPlus[jStat] = (TH1D *) gDirectory->Get(buffer);
1230  }
1231  for (int jStat=0;jStat<11;jStat++) {
1232  sprintf(buffer,"%sPMinus%s_%i", mpreFix, mKey, jStat);
1233  hPMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
1234  }
1235  for (int jStat=0;jStat<17;jStat++) {
1236  sprintf(buffer,"%sPPlusMinus%s_%i", mpreFix, mKey, jStat);
1237  hPPlusMinus[jStat] = (TH1D *) gDirectory->Get(buffer);
1238  }
1239 
1240 
1241  int mBin, oBin;
1242 
1243  for (int iPhi=1;iPhi<=mNPhiBins;iPhi++) {
1244  for (int iEta=1;iEta<=mNEtaBins;iEta++) {
1245  double NS[16], ND[2], NP[5], NM[5], NC[8];
1246  double PS[21], PD[16], PP[11], PM[11], PC[17];
1247 
1248  double sigSq0, sigSq1, sigSq2, pHat[3];
1249  double Ss[] = {0, 0, 0}, nSs = 0;
1250  double Sd[] = {0, 0, 0}, nSd = 0;
1251  double Sp[] = {0, 0, 0}, nSp = 0;
1252  double Sm[] = {0, 0, 0}, nSm = 0;
1253  double Sc1[] = {0, 0, 0}, nSc1 = 0;
1254  double Sc2[] = {0, 0, 0}, nSc2 = 0;
1255 
1256  double totEvents, sumEvents, plusEvents, minusEvents, plusMinusEvents;
1257  double nErr, pErr;
1258 
1259  mBin = (int) hnBins->GetBinContent(iPhi,iEta);
1260  oBin = (int) hoffset->GetBinContent(iPhi,iEta);
1261 
1262  int iEtaBin = 0, iPhiBin = 1;
1263  int nEtaBin = getNumEtaBins(iEta);
1264  float phi1, phi2, eta1, eta2;
1265  for (int iBin=oBin+1;iBin<=oBin+mBin;iBin++) {
1266  iEtaBin++;
1267  if (iEtaBin > nEtaBin) {
1268  iEtaBin = 1;
1269  iPhiBin++;
1270  }
1271  phi1 = 0 + 2*3.1415926*getPhiStart(iPhiBin,iPhi)/mNPhiBins;
1272  phi2 = phi1 + iPhi * 2*3.1415926 / mNPhiBins;
1273  eta1 = mEtaMin + (mEtaMax-mEtaMin)*getEtaStart(iEtaBin,iEta)/mNEtaBins;
1274  eta2 = eta1 + iEta * (mEtaMax-mEtaMin) / mNEtaBins;
1275  totEvents = hTotEvents[0]->GetBinContent(iBin);
1276  sumEvents = hTotEvents[1]->GetBinContent(iBin);
1277  if (sumEvents > 0) {
1278  for (int jStat=0;jStat<16;jStat++) {
1279  NS[jStat] = hNSum[jStat]->GetBinContent(iBin) / sumEvents;
1280  }
1281  for (int jStat=0;jStat<2;jStat++) {
1282  ND[jStat] = hNDiff[jStat]->GetBinContent(iBin) / sumEvents;
1283  }
1284  for (int jStat=0;jStat<21;jStat++) {
1285  PS[jStat] = hPSum[jStat]->GetBinContent(iBin) / sumEvents;
1286  }
1287  for (int jStat=0;jStat<16;jStat++) {
1288  PD[jStat] = hPDiff[jStat]->GetBinContent(iBin) / sumEvents;
1289  }
1290  double NSum = NS[0] + NS[1];
1291  if (NSum > 0) {
1292  if (NS[0] > 0) {
1293  pHat[1] = PS[1] / NS[0];
1294  } else {
1295  pHat[1] = 0;
1296  }
1297  if (NS[1] > 0) {
1298  pHat[2] = PS[2] / NS[1];
1299  } else {
1300  pHat[2] = 0;
1301  }
1302  sigSq0 = (PS[15] - NSum*PS[14]
1303  - pHat[1]*(NS[12] - NSum*NS[11])
1304  - pHat[2]*(NS[14] - NSum*NS[13])) / sqrt(NSum);
1305  Ss[0] += sigSq0;
1306  sigSq1 = (PS[3]+PS[4]
1307  - pHat[1]*(NS[2] + NS[3])
1308  - pHat[2]*(NS[3] + NS[4])) / NSum;
1309  Ss[1] += sigSq1;
1310  sigSq2 = (-PS[9] - PS[10]
1311  + pHat[1]*(NS[8] + NS[9])
1312  + pHat[2]*(NS[9] + NS[10])) * NSum;
1313  Ss[2] += sigSq2;
1314  nSs +=1;
1315  bTuple.type = 11;
1316  bTuple.phiScale = iPhi;
1317  bTuple.etaScale = iEta;
1318  bTuple.phi = (phi1+phi2)/2;
1319  bTuple.eta = (eta1+eta2)/2;
1320  bTuple.sig2 = sigSq0;
1321  bTuple.sig2_1 = sigSq1;
1322  bTuple.sig2_2 = sigSq2;
1323  bTuple.nbar = NSum;
1324  bTuple.events = sumEvents;
1325  bTuple.f3 = 0;
1326  binTuple->Fill(&bTuple.type);
1327 
1328  double NDiff = NS[0] - NS[1];
1329  sigSq0 = PD[10] - NDiff*PD[9]
1330  - pHat[1]*(ND[0]-NDiff*NS[11])
1331  + pHat[2]*(ND[1]-NDiff*NS[13]);
1332  Sd[0] += sigSq0 / sqrt(NSum);
1333  sigSq1 = PD[0] - PD[1]
1334  - pHat[1]*(NS[2] - NS[3])
1335  + pHat[2]*(NS[3] - NS[4]);
1336  Sd[1] += sigSq1 / NSum;
1337  sigSq2 = PD[3] - PD[4] - NDiff*(PD[6]+PD[7])
1338  - pHat[1]*(NS[5]-NS[6] - NDiff*(NS[8]+NS[9]))
1339  + pHat[2]*(NS[6]-NS[7] - NDiff*(NS[9]+NS[10]));
1340  Sd[2] += sigSq2;
1341  nSd += 1;
1342  bTuple.type = 12;
1343  bTuple.phiScale = iPhi;
1344  bTuple.etaScale = iEta;
1345  bTuple.phi = (phi1+phi2)/2;
1346  bTuple.eta = (eta1+eta2)/2;
1347  bTuple.sig2 = sigSq0/sqrt(NSum);
1348  bTuple.sig2_1 = sigSq1/NSum;
1349  bTuple.sig2_2 = sigSq2;
1350  bTuple.nbar = NSum;
1351  bTuple.events = sumEvents;
1352  bTuple.f3 = 0;
1353  binTuple->Fill(&bTuple.type);
1354  }
1355  }
1356 
1357  plusEvents = hTotEvents[2]->GetBinContent(iBin);
1358  pHat[1] = 0;
1359  if (plusEvents > 0) {
1360  for (int jStat=0;jStat<5;jStat++) {
1361  NP[jStat] = hNPlus[jStat]->GetBinContent(iBin) / plusEvents;
1362  }
1363  for (int jStat=0;jStat<11;jStat++) {
1364  PP[jStat] = hPPlus[jStat]->GetBinContent(iBin) / plusEvents;
1365  }
1366  if (NP[0] > 0) {
1367  pHat[1] = PP[1] / NP[0];
1368  sigSq0 = PP[10] - NP[0]*PP[9]
1369  - pHat[1]*NP[3] + PP[1]*NP[2];
1370  Sp[0] += sigSq0 / sqrt(NP[0]);
1371  sigSq1 = PP[2] - pHat[1]*NP[1];
1372  Sp[1] += sigSq1 / NP[0];
1373  sigSq2 = PP[1] - NP[0]*PP[3];
1374  Sp[2] += sigSq2;
1375  nSp += 1;
1376  bTuple.type = 13;
1377  bTuple.phiScale = iPhi;
1378  bTuple.etaScale = iEta;
1379  bTuple.phi = (phi1+phi2)/2;
1380  bTuple.eta = (eta1+eta2)/2;
1381  bTuple.sig2 = sigSq0/sqrt(NP[0]);
1382  bTuple.sig2_1 = sigSq1/NP[0];
1383  bTuple.sig2_2 = sigSq2;
1384  bTuple.nbar = NP[0];
1385  bTuple.events = plusEvents;
1386  bTuple.f3 = 0;
1387  binTuple->Fill(&bTuple.type);
1388  }
1389  } else {
1390  for (int jStat=0;jStat<5;jStat++) {
1391  NP[jStat] = 0;
1392  }
1393  for (int jStat=0;jStat<11;jStat++) {
1394  PP[jStat] = 0;
1395  }
1396  }
1397 
1398  minusEvents = hTotEvents[3]->GetBinContent(iBin);
1399  pHat[2] = 0;
1400  if (minusEvents > 0) {
1401  for (int jStat=0;jStat<5;jStat++) {
1402  NM[jStat] = hNMinus[jStat]->GetBinContent(iBin) / minusEvents;
1403  }
1404  for (int jStat=0;jStat<11;jStat++) {
1405  PM[jStat] = hPMinus[jStat]->GetBinContent(iBin) / minusEvents;
1406  }
1407  if (NM[0] > 0) {
1408  pHat[2] = PM[1] / NM[0];
1409  sigSq0 = PM[10] - NM[0]*PM[9]
1410  - pHat[2]*NM[3] + PM[1]*NM[2];
1411  Sm[0] += sigSq0 / sqrt(NM[0]);
1412  sigSq1 = PM[2] - pHat[2]*NM[1];
1413  Sm[1] += sigSq1 / NM[0];
1414  sigSq2 = PM[1] - NM[0]*PM[3];
1415  Sm[2] += sigSq2;
1416  nSm += 1;
1417  bTuple.type = 14;
1418  bTuple.phiScale = iPhi;
1419  bTuple.etaScale = iEta;
1420  bTuple.phi = (phi1+phi2)/2;
1421  bTuple.eta = (eta1+eta2)/2;
1422  bTuple.sig2 = sigSq0/sqrt(NM[0]);
1423  bTuple.sig2_1 = sigSq1/NM[0];
1424  bTuple.sig2_2 = sigSq2;
1425  bTuple.nbar = NM[0];
1426  bTuple.events = minusEvents;
1427  bTuple.f3 = 0;
1428  binTuple->Fill(&bTuple.type);
1429  }
1430  } else {
1431  for (int jStat=0;jStat<5;jStat++) {
1432  NM[jStat] = 0;
1433  }
1434  for (int jStat=0;jStat<11;jStat++) {
1435  PM[jStat] = 0;
1436  }
1437  }
1438 
1439 // Accumulation of hNPlusMinus required (n+ > 0) && (n- > 0)
1440 // I think this is wrong. Instead if I require (n+ > 0) || (n- > 0)
1441 // I can retrieve appropriate terms from NP and NM, at least for the
1442 // phiPt case.
1443  plusMinusEvents = hTotEvents[1]->GetBinContent(iBin);
1444  if (plusMinusEvents > 0) {
1445  for (int jStat=0;jStat<8;jStat++) {
1446  NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1447  }
1448  for (int jStat=0;jStat<17;jStat++) {
1449  PC[jStat] = hPPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1450  }
1451  sigSq0 = 0;
1452  if (NM[0] > 0) {
1453  sigSq0 = (PC[12] - NC[1]*PC[6]
1454  - pHat[1]*NC[6] + pHat[1]*NC[1]*NC[3]) / sqrt(NM[0]);
1455  Sc1[0] += sigSq0;
1456  }
1457  sigSq1 = 0;
1458  if (NP[0]*NM[0] > 0) {
1459  sigSq1 = (PC[5] - pHat[1]*NC[2]) / sqrt(NP[0]*NM[0]);
1460  Sc1[1] += sigSq1;
1461  }
1462  sigSq2 = 0;
1463  if (NM[0] > 0) {
1464  sigSq2 = (PC[8] - NC[1]*PC[2]) * sqrt(NP[0]/NM[0]);
1465  Sc1[2] += sigSq2;
1466  }
1467  nSc1 += 1;
1468  bTuple.type = 15;
1469  bTuple.phiScale = iPhi;
1470  bTuple.etaScale = iEta;
1471  bTuple.phi = (phi1+phi2)/2;
1472  bTuple.eta = (eta1+eta2)/2;
1473  bTuple.sig2 = sigSq0;
1474  bTuple.sig2_1 = sigSq1;
1475  bTuple.sig2_2 = sigSq2;
1476  bTuple.nbar = sqrt(NM[0]*NP[0]);
1477  bTuple.events = plusMinusEvents;
1478  bTuple.f3 = 0;
1479  binTuple->Fill(&bTuple.type);
1480 
1481  sigSq0 = 0;
1482  if (NP[0] > 0) {
1483  sigSq0 = (PC[13] - NC[0]*PC[7]
1484  - pHat[2]*NC[7] + pHat[2]*NC[0]*NC[4]) / sqrt(NP[0]);
1485  Sc2[0] += sigSq0;
1486  }
1487  sigSq1 = 0;
1488  if (NP[0]*NM[0] > 0) {
1489  sigSq1 = (PC[4] - pHat[2]*NC[2]) / sqrt(NP[0]*NM[0]);
1490  Sc2[1] += sigSq1;
1491  }
1492  sigSq2 = 0;
1493  if (NP[0] > 0) {
1494  sigSq2 = (PC[9] - NC[0]*PC[3]) * sqrt(NM[0]/NP[0]);
1495  Sc2[2] += sigSq2;
1496  }
1497  nSc2 += 1;
1498  bTuple.type = 16;
1499  bTuple.phiScale = iPhi;
1500  bTuple.etaScale = iEta;
1501  bTuple.phi = (phi1+phi2)/2;
1502  bTuple.eta = (eta1+eta2)/2;
1503  bTuple.sig2 = sigSq0;
1504  bTuple.sig2_1 = sigSq1;
1505  bTuple.sig2_2 = sigSq2;
1506  bTuple.nbar = sqrt(NM[0]*NP[0]);
1507  bTuple.events = plusMinusEvents;
1508  bTuple.f3 = 0;
1509  binTuple->Fill(&bTuple.type);
1510  }
1511 // plusMinusEvents = hTotEvents[4]->GetBinContent(iBin);
1512 // if (plusMinusEvents > 0) {
1513 // for (int jStat=0;jStat<8;jStat++) {
1514 // NC[jStat] = hNPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1515 // }
1516 // for (int jStat=0;jStat<17;jStat++) {
1517 // PC[jStat] = hPPlusMinus[jStat]->GetBinContent(iBin) / plusMinusEvents;
1518 // }
1519 // if (NC[0] > 0) {
1520 // pHat[1] = PC[0] / NC[0];
1521 // } else {
1522 // pHat[1] = 0;
1523 // }
1524 // if (NC[1] > 0) {
1525 // pHat[2] = PC[1] / NC[1];
1526 // } else {
1527 // pHat[2] = 0;
1528 // }
1529 //
1530 // sigSq = (PC[12] - NC[1]*PC[6]
1531 // - pHat[1]*NC[6] + pHat[1]*NC[1]*NC[3]) / sqrt(NC[1]);
1532 // Sc1[0] += sigSq;
1533 // sigSq = (PC[5] - NC[1]*PC[0]
1534 // - pHat[1]*NC[2] + PC[0]*NC[1]) / sqrt(NC[0]*NC[1]);
1535 // Sc1[1] += sigSq;
1536 // sigSq = (PC[8] - NC[1]*PC[2]) * sqrt(NC[0]/NC[1]);
1537 // Sc1[2] += sigSq;
1538 // nSc1 += 1;
1539 //
1540 // sigSq = (PC[13] - NC[0]*PC[7]
1541 // - pHat[2]*NC[7] + pHat[2]*NC[0]*NC[4]) / sqrt(NC[0]);
1542 // Sc2[0] += sigSq;
1543 // sigSq = (PC[4] - NC[0]*PC[1]
1544 // - pHat[2]*NC[2] + PC[1]*NC[0]) / sqrt(NC[0]*NC[1]);
1545 // Sc2[1] += sigSq;
1546 // sigSq = (PC[9] - NC[0]*PC[3]) * sqrt(NC[1]/NC[0]);
1547 // Sc2[2] += sigSq;
1548 // nSc2 += 1;
1549 // }
1550  }
1551 
1552  if (nSs > 0) {
1553  PNSig[0]->SetBinContent(iPhi,iEta,Ss[0]/nSs);
1554  PNSig[1]->SetBinContent(iPhi,iEta,Ss[1]/nSs);
1555  PNSig[2]->SetBinContent(iPhi,iEta,Ss[2]/nSs);
1556  }
1557  if (nSd > 0) {
1558  PNDel[0]->SetBinContent(iPhi,iEta,Sd[0]/nSd);
1559  PNDel[1]->SetBinContent(iPhi,iEta,Sd[1]/nSd);
1560  PNDel[2]->SetBinContent(iPhi,iEta,Sd[2]/nSd);
1561  }
1562  if (nSp > 0) {
1563  PNPlus[0]->SetBinContent(iPhi,iEta,Sp[0]/nSp);
1564  PNPlus[1]->SetBinContent(iPhi,iEta,Sp[1]/nSp);
1565  PNPlus[2]->SetBinContent(iPhi,iEta,Sp[2]/nSp);
1566  }
1567  if (nSm > 0) {
1568  PNMinus[0]->SetBinContent(iPhi,iEta,Sm[0]/nSm);
1569  PNMinus[1]->SetBinContent(iPhi,iEta,Sm[1]/nSm);
1570  PNMinus[2]->SetBinContent(iPhi,iEta,Sm[2]/nSm);
1571  }
1572  if (nSc1) {
1573  PNPlusMinus[0]->SetBinContent(iPhi,iEta,Sc1[0]/nSc1);
1574  PNPlusMinus[1]->SetBinContent(iPhi,iEta,Sc1[1]/nSc1);
1575  PNPlusMinus[2]->SetBinContent(iPhi,iEta,Sc1[2]/nSc1);
1576  }
1577  if (nSc2) {
1578  PNMinusPlus[0]->SetBinContent(iPhi,iEta,Sc2[0]/nSc2);
1579  PNMinusPlus[1]->SetBinContent(iPhi,iEta,Sc2[1]/nSc2);
1580  PNMinusPlus[2]->SetBinContent(iPhi,iEta,Sc2[2]/nSc2);
1581  }
1582 // I haven't been able to find a reasonable approximation to the errors
1583 // on the Pt-N covariance terms. My best guess is the errors factor,
1584 // and since we already have the N and the Pt errors I use those.
1585  nErr = NSigErrors->GetBinContent(iPhi,iEta);
1586  pErr = PSigErrors->GetBinContent(iPhi,iEta);
1587  PNSigErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1588  nErr = NDelErrors->GetBinContent(iPhi,iEta);
1589  pErr = PDelErrors->GetBinContent(iPhi,iEta);
1590  PNDelErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1591  nErr = NPlusErrors->GetBinContent(iPhi,iEta);
1592  pErr = PPlusErrors->GetBinContent(iPhi,iEta);
1593  PNPlusErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1594  nErr = NMinusErrors->GetBinContent(iPhi,iEta);
1595  pErr = PMinusErrors->GetBinContent(iPhi,iEta);
1596  PNMinusErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1597  nErr = NPlusMinusErrors->GetBinContent(iPhi,iEta);
1598  pErr = PPlusMinusErrors->GetBinContent(iPhi,iEta);
1599  PNPlusMinusErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1600  PNMinusPlusErrors->SetBinContent(iPhi,iEta,sqrt(nErr*pErr));
1601  }
1602  }
1603 }
1604 
1605 //
1606 //------------ Below are init, delete, write functions -------///
1607 //
1608 
1609 
1610 //--------------------------------------------------------------------------
1611 void StEStructSigmas::writeHistograms() {
1612  NSig->Write();
1613  NDel->Write();
1614  NPlus->Write();
1615  NMinus->Write();
1616  NPlusMinus->Write();
1617  NSigCorrection->Write();
1618  NDelCorrection->Write();
1619  NPlusCorrection->Write();
1620  NMinusCorrection->Write();
1621  NPlusMinusCorrection->Write();
1622  for (int iType=0;iType<3;iType++) {
1623  PSig[iType]->Write();
1624  PDel[iType]->Write();
1625  PPlus[iType]->Write();
1626  PMinus[iType]->Write();
1627  PPlusMinus[iType]->Write();
1628  PNSig[iType]->Write();
1629  PNDel[iType]->Write();
1630  PNPlus[iType]->Write();
1631  PNMinus[iType]->Write();
1632  PNPlusMinus[iType]->Write();
1633  PNMinusPlus[iType]->Write();
1634  }
1635  PSig[3]->Write();
1636 
1637  SPtHat->Write();
1638  PPtHat->Write();
1639  MPtHat->Write();
1640  sigSPtHat->Write();
1641  sigPPtHat->Write();
1642  sigMPtHat->Write();
1643 
1644  NSigErrors->Write();
1645  NDelErrors->Write();
1646  NPlusErrors->Write();
1647  NMinusErrors->Write();
1648  NPlusMinusErrors->Write();
1649  PSigErrors->Write();
1650  PDelErrors->Write();
1651  PPlusErrors->Write();
1652  PMinusErrors->Write();
1653  PPlusMinusErrors->Write();
1654  PNSigErrors->Write();
1655  PNDelErrors->Write();
1656  PNPlusErrors->Write();
1657  PNMinusErrors->Write();
1658  PNPlusMinusErrors->Write();
1659  PNMinusPlusErrors->Write();
1660 
1661  sigSPtHatErrors->Write();
1662  sigPPtHatErrors->Write();
1663  sigMPtHatErrors->Write();
1664  binTuple->Write();
1665  scaleTuple->Write();
1666  sumTuple->Write();
1667 }
1668 
1669 //--------------------------------------------------------------------------
1670 void StEStructSigmas::initHistograms() {
1671  char line[255];
1672 
1673  // Here are histograms for storing variances and errors of variances..
1674  sprintf( line, "NSig_%s", mKey );
1675  NSig = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1676  sprintf( line, "NDel_%s", mKey );
1677  NDel = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1678  sprintf( line, "NPlus_%s", mKey );
1679  NPlus = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1680  sprintf( line, "NMinus_%s", mKey );
1681  NMinus = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1682  sprintf( line, "NPlusMinus_%s", mKey );
1683  NPlusMinus = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1684  sprintf( line, "NSigCorrection_%s", mKey );
1685  NSigCorrection = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1686  sprintf( line, "NDelCorrection_%s", mKey );
1687  NDelCorrection = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1688  sprintf( line, "NPlusCorrection_%s", mKey );
1689  NPlusCorrection = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1690  sprintf( line, "NMinusCorrection_%s", mKey );
1691  NMinusCorrection = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1692  sprintf( line, "NPlusMinusCorrection_%s", mKey );
1693  NPlusMinusCorrection = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1694 
1695  for (int iType=0;iType<3;iType++) {
1696  sprintf( line, "PSig%i_%s", iType, mKey );
1697  PSig[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1698  sprintf( line, "PDel%i_%s", iType, mKey );
1699  PDel[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1700  sprintf( line, "PPlus%i_%s", iType, mKey );
1701  PPlus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1702  sprintf( line, "PMinus%i_%s", iType, mKey );
1703  PMinus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1704  sprintf( line, "PPlusMinus%i_%s", iType, mKey );
1705  PPlusMinus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1706  sprintf( line, "PNSig%i_%s", iType, mKey );
1707  PNSig[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1708  sprintf( line, "PNDel%i_%s", iType, mKey );
1709  PNDel[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1710  sprintf( line, "PNPlus%i_%s", iType, mKey );
1711  PNPlus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1712  sprintf( line, "PNMinus%i_%s", iType, mKey );
1713  PNMinus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1714  sprintf( line, "PNPlusMinus%i_%s", iType, mKey );
1715  PNPlusMinus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1716  sprintf( line, "PNMinusPlus%i_%s", iType, mKey );
1717  PNMinusPlus[iType] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1718  }
1719  sprintf( line, "PSig%i_%s", 3, mKey );
1720  PSig[3] = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1721 
1722  sprintf( line, "SptHat_%s", mKey );
1723  SPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1724  sprintf( line, "PptHat_%s", mKey );
1725  PPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1726  sprintf( line, "MptHat_%s", mKey );
1727  MPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1728  sprintf( line, "SsigPtHat_%s", mKey );
1729  sigSPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1730  sprintf( line, "PsigPtHat_%s", mKey );
1731  sigPPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1732  sprintf( line, "MsigPtHat_%s", mKey );
1733  sigMPtHat = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1734 
1735  sprintf( line, "NSigErrors_%s", mKey );
1736  NSigErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1737  sprintf( line, "NDelErrors_%s", mKey );
1738  NDelErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1739  sprintf( line, "NPlusErrors_%s", mKey );
1740  NPlusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1741  sprintf( line, "NMinusErrors_%s", mKey );
1742  NMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1743  sprintf( line, "NPlusMinusErrors_%s", mKey );
1744  NPlusMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1745  sprintf( line, "PSigErrors_%s", mKey );
1746  PSigErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1747  sprintf( line, "PDelErrors_%s", mKey );
1748  PDelErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1749  sprintf( line, "PPlusErrors_%s", mKey );
1750  PPlusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1751  sprintf( line, "PMinusErrors_%s", mKey );
1752  PMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1753  sprintf( line, "PPlusMinusErrors_%s", mKey );
1754  PPlusMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1755  sprintf( line, "PNSigErrors_%s", mKey );
1756  PNSigErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1757  sprintf( line, "PNDelErrors_%s", mKey );
1758  PNDelErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1759  sprintf( line, "PNPlusErrors_%s", mKey );
1760  PNPlusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1761  sprintf( line, "PNMinusErrors_%s", mKey );
1762  PNMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1763  sprintf( line, "PNPlusMinusErrors_%s", mKey );
1764  PNPlusMinusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1765  sprintf( line, "PNMinusPlusErrors_%s", mKey );
1766  PNMinusPlusErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1767 
1768  sprintf( line, "SsigPtHatErrors_%s", mKey );
1769  sigSPtHatErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1770  sprintf( line, "PsigPtHatErrors_%s", mKey );
1771  sigPPtHatErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1772  sprintf( line, "MsigPtHatErrors_%s", mKey );
1773  sigMPtHatErrors = new TH2D(line,line,mNPhiBins,0.0,6.2831852,mNEtaBins,0.0,mEtaMax-mEtaMin);
1774 
1775  sprintf( line, "binTuple_%s", mKey );
1776  binTuple = new TNtuple(line,"sigma dependence on eta",
1777  "type:phiScale:etaScale:phi:eta:sig2:sig2_1:sig2_2:nbar:events:f3");
1778  sprintf( line, "scaleTuple_%s", mKey );
1779  scaleTuple = new TNtuple(line,"fitting dependence on eta",
1780  "type:phiScale:etaScale:A:B:nBins:f3:f3sq:sig2:sig2f3");
1781  sprintf( line, "sumTuple_%s", mKey );
1782  sumTuple = new TNtuple(line,"fit parameters for eta",
1783  "type:B:nBins");
1784 }
1785 
1786 //--------------------------------------------------------------------------
1787 void StEStructSigmas::deleteHistograms() {
1788 
1789 printf("Deleting bin***Var2, bin***Errors histograms.\n");
1790  delete NSig;
1791  delete NDel;
1792  delete NPlus;
1793  delete NMinus;
1794  delete NPlusMinus;
1795  delete NSigCorrection;
1796  delete NDelCorrection;
1797  delete NPlusCorrection;
1798  delete NMinusCorrection;
1799  delete NPlusMinusCorrection;
1800  for (int iType=0;iType<3;iType++) {
1801  delete PSig[iType];
1802  delete PDel[iType];
1803  delete PPlus[iType];
1804  delete PMinus[iType];
1805  delete PPlusMinus[iType];
1806  delete PNSig[iType];
1807  delete PNDel[iType];
1808  delete PNPlus[iType];
1809  delete PNMinus[iType];
1810  delete PNPlusMinus[iType];
1811  delete PNMinusPlus[iType];
1812  }
1813  delete PSig[3];
1814 
1815  delete SPtHat;
1816  delete PPtHat;
1817  delete MPtHat;
1818  delete sigSPtHat;
1819  delete sigPPtHat;
1820  delete sigMPtHat;
1821 
1822  delete NSigErrors;
1823  delete NDelErrors;
1824  delete NPlusErrors;
1825  delete NMinusErrors;
1826  delete NPlusMinusErrors;
1827  delete PSigErrors;
1828  delete PDelErrors;
1829  delete PPlusErrors;
1830  delete PMinusErrors;
1831  delete PPlusMinusErrors;
1832  delete PNSigErrors;
1833  delete PNDelErrors;
1834  delete PNPlusErrors;
1835  delete PNMinusErrors;
1836  delete PNPlusMinusErrors;
1837  delete PNMinusPlusErrors;
1838 
1839  delete sigSPtHatErrors;
1840  delete sigPPtHatErrors;
1841  delete sigMPtHatErrors;
1842 
1843  delete binTuple;
1844  delete scaleTuple;
1845  delete sumTuple;
1846 }
1847 //
1848 // Copy these from StEStructFluctAnal.cxx.
1849 // Probably should declare these routines static
1850 // and leave them in that file.
1851 //
1852 // iEta runs from 1 up to the number of etaBins that fit in,
1853 // according to the binning mode.
1854 // Return starting bin (first bin is called 0.)
1855 // When iEta is too big return -1.
1856 int StEStructSigmas::getEtaStart( int iEta, int dEta ) {
1857  if (dEta > NETABINS) {
1858  return -1;
1859  }
1860  if (1 == mEtaSumMode) {
1861  int nEta = NETABINS / dEta;
1862  int oEta = NETABINS % dEta;
1863  if ( iEta*dEta <= NETABINS ) {
1864  return (iEta-1)*dEta;
1865  }
1866  if (0 == oEta) {
1867  return -1;
1868  }
1869  if (oEta+(iEta-nEta)*dEta <= NETABINS) {
1870  return oEta+(iEta-nEta-1)*dEta;
1871  }
1872  return -1;
1873  } else if (2 == mEtaSumMode) {
1874  if (iEta+dEta <= NETABINS) {
1875  return iEta - 1;
1876  } else {
1877  return -1;
1878  }
1879  } else if (3 == mEtaSumMode) {
1880  if (iEta > 1) {
1881  return -1;
1882  } else {
1883  return 0;
1884  }
1885  } else if (4 == mEtaSumMode) {
1886  if (iEta > 1) {
1887  return -1;
1888  } else {
1889  return NETABINS-dEta-1;
1890  }
1891  } else if (5 == mEtaSumMode) {
1892  if (iEta > 1) {
1893  return -1;
1894  } else {
1895  return (NETABINS-dEta) / 2;
1896  }
1897  } else if (6 == mEtaSumMode) {
1898  int nEta = NETABINS / dEta;
1899  if (iEta > nEta) {
1900  return -1;
1901  }
1902  int oEta = (NETABINS % dEta) / 2;
1903  return oEta + (iEta-1)*dEta;
1904  }
1905  return -1;
1906 }
1907 int StEStructSigmas::getPhiStart( int iPhi, int dPhi ) {
1908  if (dPhi > NPHIBINS) {
1909  return -1;
1910  }
1911  if (1 == mPhiSumMode) {
1912  int nPhi = NPHIBINS / dPhi;
1913  int oPhi = NPHIBINS % dPhi;
1914  if ( iPhi*dPhi <= NPHIBINS ) {
1915  return (iPhi-1)*dPhi;
1916  }
1917  if (0 == oPhi) {
1918  return -1;
1919  }
1920  if (oPhi+(iPhi-nPhi)*dPhi <= NPHIBINS) {
1921  return oPhi+(iPhi-nPhi-1)*dPhi;
1922  }
1923  return -1;
1924  } else if (2 == mPhiSumMode) {
1925  if (iPhi+dPhi <= NPHIBINS) {
1926  return iPhi - 1;
1927  } else {
1928  return -1;
1929  }
1930  } else if (3 == mPhiSumMode) {
1931  if (iPhi > 1) {
1932  return -1;
1933  } else {
1934  return 0;
1935  }
1936  } else if (4 == mPhiSumMode) {
1937  if (iPhi > 1) {
1938  return -1;
1939  } else {
1940  return NPHIBINS-dPhi - 1;
1941  }
1942  } else if (5 == mPhiSumMode) {
1943  if (iPhi > 1) {
1944  return -1;
1945  } else {
1946  return (NPHIBINS-dPhi) / 2;
1947  }
1948  } else if (6 == mPhiSumMode) {
1949  int nPhi = NPHIBINS / dPhi;
1950  if (iPhi > nPhi) {
1951  return -1;
1952  }
1953  int oPhi = (NPHIBINS % dPhi) / 2;
1954  return oPhi + (iPhi-1)*dPhi;
1955  }
1956  return -1;
1957 }
1958 int StEStructSigmas::getNumEtaBins( int dEta ) {
1959  if (1 == mEtaSumMode) {
1960  int nEta = NETABINS / dEta;
1961  int oEta = NETABINS % dEta;
1962  if ( 0 == oEta ) {
1963  return nEta;
1964  } else {
1965  return 2 * nEta;
1966  }
1967  } else if (2 == mEtaSumMode) {
1968  return NETABINS + 1 - dEta;
1969  } else if (3 == mEtaSumMode) {
1970  return 1;
1971  } else if (4 == mEtaSumMode) {
1972  return 1;
1973  } else if (5 == mEtaSumMode) {
1974  return 1;
1975  } else if (5 == mEtaSumMode) {
1976  int nEta = NETABINS / dEta;
1977  return nEta;
1978  }
1979  return 0;
1980 }
1981 int StEStructSigmas::getNumPhiBins( int dPhi ) {
1982  if (1 == mPhiSumMode) {
1983  int nPhi = NPHIBINS / dPhi;
1984  int oPhi = NPHIBINS % dPhi;
1985  if ( 0 == oPhi ) {
1986  return nPhi;
1987  } else {
1988  return 2 * nPhi;
1989  }
1990  } else if (2 == mPhiSumMode) {
1991  return NPHIBINS + 1 - dPhi;
1992  } else if (3 == mPhiSumMode) {
1993  return 1;
1994  } else if (4 == mPhiSumMode) {
1995  return 1;
1996  } else if (5 == mPhiSumMode) {
1997  return 1;
1998  } else if (5 == mPhiSumMode) {
1999  int nPhi = NPHIBINS / dPhi;
2000  return nPhi;
2001  }
2002  return 0;
2003 }
Generated by   doxygen 1.8.5