StRoot  1
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
StEStructHAdd.cxx
1 /**********************************************************************
2  *
3  * $Id: StEStructHAdd.cxx,v 1.17 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 "StEStructHAdd.h"
15 
16 #include "StEStructPool/Correlations/StEStructBinning.h"
17 #include "StEStructPool/Correlations/StEStructCutBin.h"
18 
19 #include "Stiostream.h"
20 #include "TH1.h"
21 #include "TH2.h"
22 #include "TH2D.h"
23 #include "TFile.h"
24 
25 ClassImp(StEStructHAdd)
26 
27  /*
28  * 12/03/08 djp Instrumented StEStruct2ptCorrelations so we can select specific
29  * groups of histograms. We still try to read all the possible histograms
30  * but don't get upset if we don't find them.
31  * Expunge code to look for old style of name (without zBin) to make code simpler.
32  */
33  /*
34  * 06/05/08 djp Looking at SEtaDPhi histograms I was confused that they were not symmetric
35  * around SEta=0. I guess this was intentional with the idea that forward and
36  * backward could be different. Should re-think this as it may make more
37  * sense to symmetrize.
38  */
39  /* 09/27/12 djp Now taking eta range from Cuts file. Had been interrogating
40  * binning class for detaMax. Extract from special purpose histogram instead.
41  */
42  //------------------------------------------------------------------------
43 void StEStructHAdd::addCuts(const char* outfile, TFile* inFile,
44  int* nlist, int ntot, int parentDist[][2], int* nParentDist, int symmXX) {
45 
46  const char* base[]={"Sib","Mix"};
47  const char* tpe[]={"pp","pm","mp","mm"};
48  const char* knd[]={"YtYt", "NYtYt", "PtPt",
49  "PhiPhi", "NPhiPhi", "PrPhiPhi", "PaPhiPhi", "PbPhiPhi",
50  "EtaEta", "PrEtaEta", "PaEtaEta", "PbEtaEta",
51  "EtaEtaSS", "PrEtaEtaSS", "PaEtaEtaSS", "PbEtaEtaSS",
52  "EtaEtaAS", "PrEtaEtaAS", "PaEtaEtaAS", "PbEtaEtaAS",
53  "DEtaDPhiArr", "NDEtaDPhiArr", "PrDEtaDPhiArr", "PaDEtaDPhiArr", "PbDEtaDPhiArr",
54  "SYtDYt", "NSYtDYt",
55  "SEtaDPhiArr", "NSEtaDPhiArr", "PrSEtaDPhiArr", "PaSEtaDPhiArr", "PbSEtaDPhiArr",
56  "Qinv", "NQinv"};
57  int isXXHist[]={1, 1, 1,
58  1, 1, 1, 1, 1,
59  1, 1, 1, 1,
60  1, 1, 1, 1,
61  1, 1, 1, 1,
62  0, 0, 0, 0, 0,
63  0, 0,
64  0, 0, 0, 0, 0,
65  0, 0};
66  const char* symEtaPhi[]={"DEtaDPhi", "NDEtaDPhi", "PrDEtaDPhi", "PaDEtaDPhi", "PbDEtaDPhi"};
67  const char* symPhi[]={"SEtaDPhi", "NSEtaDPhi", "PrSEtaDPhi", "PaSEtaDPhi", "PbSEtaDPhi"};
68  const char* Title[]={"Sibling","Mixed"};
69  const char* Species[]={"A","B"};
70  const char* Type[]={" : +.+"," : +.-"," : -.+"," : -.-"};
71  StEStructBinning* b=StEStructBinning::Instance();
72  TH2D *hEtaPhi;
73  inFile->GetObject("EtaPhiRange",hEtaPhi);
74  double etaMin = -1;
75  double etaMax = +1;
76  if (hEtaPhi) {
77  TAxis *x = hEtaPhi->GetXaxis();
78  etaMin = x->GetXmin();
79  etaMax = x->GetXmax();
80  }
81  b->setEtaRange(etaMin,etaMax);
82 
83  StEStructCutBin* cb = StEStructCutBin::Instance();
84 
85  TFile* outFile = new TFile(outfile,"RECREATE");
86  outFile->cd();
87  if (hEtaPhi) {
88  hEtaPhi->Write();
89  }
90 
91  for (int i=0;i<2;i++) {
92  for (int j=0;j<4;j++) {
93  for (int k=0;k<34;k++) {
94 
95  TH2D *outhist=0;
96  TH2D *tmp=0, *cpy=0;
97  inFile->cd();
98  int zBin = 0;
99  while (zBin<99) {
100  // Histogram name is of form
101  // SibppDEtaDPhi_cutBin_n_zBuf_m
102  // where n is the usual cut bin number and m is the zbuffer index.
103  TString htype(base[i]); htype+=tpe[j]; htype+=knd[k];
104  for (int n=0;n<ntot;n++) {
105  TString fullName(htype.Data()); fullName+="_cutBin_"; fullName+=nlist[n]; fullName+="_zBuf_"; fullName+=zBin;
106  TString hName(htype.Data()); hName+="_zBuf_"; hName+=zBin;
107  inFile->GetObject(fullName.Data(),tmp);
108  if (!tmp) {
109  // Have either processed all zbins or this histogram (at least for zBuf=0) isn't in the file.
110  goto lastZ;
111  }
112  if (0==n) {
113  outhist=(TH2D *)tmp->Clone();
114  if (symmXX && isXXHist[k]) {
115  // This part is (was) only intended to be used with mode 5 (pid)
116  if (cb->notSymmetrizedXX(nlist[n],j)) {
117  symmetrizeXX(outhist);
118  }
119  }
120  outhist->SetName(hName.Data());
121  outhist->SetTitle(tmp->GetTitle());
122  } else {
123  cpy = (TH2D *)tmp->Clone();
124  if (symmXX && isXXHist[k]) {
125  // This part is (was) only intended to be used with mode 5 (pid)
126  if (cb->notSymmetrizedXX(nlist[n],j)) {
127  symmetrizeXX(cpy);
128  }
129  }
130  outhist->Add(cpy);
131  delete cpy;
132  }
133  }
134  for (int isym=0;isym<5;isym++) {
135  if (!strncmp(symEtaPhi[isym],knd[k],8)) {
136  // Here outhist is not symmetrized. It will be replaced by tmp which is.
137  b->setNDEtaBins(outhist->GetNbinsX());
138  b->setNDPhiBins(outhist->GetNbinsY());
139  TString hname(base[i]); hname+=tpe[j]; hname+=symEtaPhi[isym];
140  TString zName(hname); zName+="_zBuf_"; zName+=zBin;
141  TString htitle(Title[i]); htitle+=Type[j]; htitle+=symEtaPhi[isym];
142  TString zTitle(htitle); zTitle+="_zBuf_"; zTitle+=zBin;
143  tmp = new TH2D(zName.Data(),zTitle.Data(),b->hdetaBins(),b->detaMin(),b->detaMax(),b->hdphiBins(),b->dphiMin(),b->dphiMax());
144  // Turn error propagation on by default since it seems setting bin content bypasses that.
145  // For a few histogram types we handle errors explicitly to avoid double counting.
146  tmp->Sumw2();
147  for(int ieta=0;ieta<b->detaBins();ieta++){
148  for(int iphi=0;iphi<b->dphiBins();iphi++){
149  float eta1 = b->detaVal(ieta);
150  float eta2 = -eta1;
151  int ieta1 = b->hdetaBin(eta1);
152  int ieta2 = b->hdetaBin(eta2);
153  float phi1 = b->dphiVal(iphi,1);
154  float phi2 = b->dphiVal(iphi,2);
155  int iphi1 = b->hdphiBin(phi1);
156  int iphi2 = b->hdphiBin(phi2);
157  double val = outhist->GetBinContent(ieta+1,iphi+1);
158  double err = outhist->GetBinError(ieta+1,iphi+1);
159  tmp->SetBinContent(ieta1,iphi1,val);
160  tmp->SetBinContent(ieta2,iphi1,val);
161  tmp->SetBinContent(ieta1,iphi2,val);
162  tmp->SetBinContent(ieta2,iphi2,val);
163  tmp->SetBinError(ieta1,iphi1,err);
164  tmp->SetBinError(ieta2,iphi1,err);
165  tmp->SetBinError(ieta1,iphi2,err);
166  tmp->SetBinError(ieta2,iphi2,err);
167  }
168  }
169  // Adjust for bins at \eta_\Delta=0 and \phi_\Delta=0 and pi
170  // Only do the adjustment if bin is centered.
171  // (When calculating \Delta\rho/\rho_{ref} this adjustment divides out,
172  // but when looking at intermediate histograms it should be less surprising.)
173  float phi1 = b->dphiVal(0,1);
174  float phi2 = b->dphiVal(0,2);
175  int iphi1 = b->hdphiBin(phi1);
176  int iphi2 = b->hdphiBin(phi2);
177  if (iphi1 == iphi2) {
178  for (int ieta=1;ieta<=b->hdetaBins();ieta++) {
179  double val = 2*tmp->GetBinContent(ieta,iphi1);
180  double err = 2*tmp->GetBinError(ieta,iphi1);
181  tmp->SetBinContent(ieta,iphi1,val);
182  tmp->SetBinError(ieta,iphi1,err);
183  }
184  }
185  // indices in StEStructBinning start at 0, histogram bins start at 1
186  phi1 = b->dphiVal(outhist->GetNbinsY()-1,1);
187  phi2 = b->dphiVal(outhist->GetNbinsY()-1,2);
188  iphi1 = b->hdphiBin(phi1);
189  iphi2 = b->hdphiBin(phi2);
190  if (iphi1 == iphi2) {
191  for (int ieta=1;ieta<=b->hdetaBins();ieta++) {
192  double val = 2*tmp->GetBinContent(ieta,iphi1);
193  double err = 2*tmp->GetBinError(ieta,iphi1);
194  tmp->SetBinContent(ieta,iphi1,val);
195  tmp->SetBinError(ieta,iphi1,err);
196  }
197  }
198  float eta1 = b->detaVal(0);
199  float eta2 = -eta1;
200  int ieta1 = b->hdetaBin(eta1);
201  int ieta2 = b->hdetaBin(eta2);
202  if (ieta1 == ieta2) {
203  for (int iphi=1;iphi<=b->hdphiBins();iphi++) {
204  double val = 2*tmp->GetBinContent(ieta1,iphi);
205  double err = 2*tmp->GetBinError(ieta1,iphi);
206  tmp->SetBinContent(ieta1,iphi,val);
207  tmp->SetBinError(ieta1,iphi,err);
208  }
209  }
210  // fill the repeated \phi_\Delta bin
211  if (0 == tmp->GetBinContent(1,1)) {
212  int iphi = b->hdphiBins();
213  for (int ieta=1;ieta<=b->hdetaBins();ieta++) {
214  double val = tmp->GetBinContent(ieta,iphi);
215  double err = tmp->GetBinError(ieta,iphi);
216  tmp->SetBinContent(ieta,1,val);
217  tmp->SetBinError(ieta,1,err);
218  }
219  }
220  // Create rotated view from -pi to pi for CD
221  hname+="Rot"; hname +="_zBuf_", hname+=zBin;
222  htitle+="Rot"; htitle+="_zBuf_", htitle+=zBin;
223  TH2D* rot = new TH2D(hname.Data(),htitle.Data(),b->hdetaBins(),b->detaMin(),b->detaMax(),b->hdphiBins(),-M_PI-M_PI/24,M_PI+M_PI/24);
224  for(int ieta=0;ieta<b->detaBins();ieta++){
225  for(int iphi=0;iphi<b->dphiBins();iphi++) {
226  float eta1 = b->detaVal(ieta);
227  float eta2 = -eta1;
228  float phi1 = b->dphiVal(iphi,1);
229  if (phi1>M_PI) phi1=2*M_PI - phi1; // get phi1 to be [0..pi]
230  float phi2 = -phi1; // phi2 [-pi..0]
231  double val = outhist->GetBinContent(ieta+1,iphi+1);
232  rot->Fill(eta1,phi1,val);
233  rot->Fill(eta2,phi1,val);
234  rot->Fill(eta1,phi2,val);
235  rot->Fill(eta2,phi2,val);
236  }
237  }
238  int bin1 = rot->GetYaxis()->FindBin(M_PI); // need to double bin values at pi and -pi
239  int bin2 = rot->GetYaxis()->FindBin(-M_PI);
240  for(int ieta=0;ieta<rot->GetNbinsX();ieta++) {
241  rot->SetBinContent(ieta+1,bin1, 2*rot->GetBinContent(ieta+1,bin1));
242  rot->SetBinContent(ieta+1,bin2, 2*rot->GetBinContent(ieta+1,bin2));
243  }
244  // Note that we need to adjust errors here also.
245  outFile->cd();
246  rot->Write();
247  delete rot;
248  // end of rotated view
249  delete outhist;
250  outhist = tmp;
251  }
252  }
253  for (int isym=0;isym<5;isym++) {
254  if (!strncmp(symPhi[isym],knd[k],8)) {
255  b->setNSEtaBins(outhist->GetNbinsX());
256  b->setNDPhiBins(outhist->GetNbinsY());
257  TString hname(base[i]);
258  hname+=tpe[j];
259  hname+=symPhi[isym];
260  TString zName(hname);
261  zName+="_zBuf_";
262  zName+=zBin;
263  TString htitle(Title[i]);
264  htitle+=Type[j];
265  htitle+=symPhi[isym];
266  TString zTitle(htitle);
267  zTitle+="_zBuf_";
268  zTitle+=zBin;
269  tmp = new TH2D(zName.Data(),zTitle.Data(),b->hdetaBins(),b->detaMin(),b->detaMax(),b->hdphiBins(),b->dphiMin(),b->dphiMax());
270  for(int ieta=0;ieta<b->setaBins();ieta++){
271  for(int iphi=0;iphi<b->dphiBins();iphi++){
272  float eta1 = b->setaVal(ieta);
273  int ieta1 = b->hdetaBin(eta1);
274  float phi1 = b->dphiVal(iphi,1);
275  float phi2 = b->dphiVal(iphi,2);
276  int iphi1 = b->hdphiBin(phi1);
277  int iphi2 = b->hdphiBin(phi2);
278  double val = outhist->GetBinContent(ieta+1,iphi+1);
279  double err = outhist->GetBinError(ieta+1,iphi+1);
280  tmp->SetBinContent(ieta1,iphi1,val);
281  tmp->SetBinContent(ieta1,iphi2,val);
282  tmp->SetBinError(ieta1,iphi1,err);
283  tmp->SetBinError(ieta1,iphi2,err);
284  }
285  }
286  if (tmp->GetBinContent(1,1)==0) { // fill the repeated bin
287  for (int ieta=0;ieta<tmp->GetNbinsX();ieta++) {
288  double val = tmp->GetBinContent(ieta+1,tmp->GetNbinsY());
289  double err = tmp->GetBinError(ieta+1,tmp->GetNbinsY());
290  tmp->SetBinContent(ieta+1,1, val);
291  tmp->SetBinError(ieta+1,1, err);
292  }
293  }
294  delete outhist;
295  outhist = tmp;
296  }
297  }
298  outFile->cd();
299  if (outhist) {
300  outhist->Write();
301  delete outhist;
302  }
303  zBin++;
304  }
305  lastZ:
306  continue;
307  }
308  }
309  }
310 
311 
312  inFile->cd();
313  TH1D* sib; inFile->GetObject("NEventsSame",sib);
314  TH1D* mix; inFile->GetObject("NEventsMixed",mix);
315  if(sib && mix) {
316  // Histogram with older naming style was found.
317  // No z-binning here.
318  outFile->cd();
319  sib->Write();
320  mix->Write();
321  } else {
322  int zBin = 0;
323  TString hSibName("NEventsSib_zBuf_");
324  TString hMixName("NEventsMix_zBuf_");
325  TString hSibPosName("NEventsPosSib_zBuf_");
326  TString hMixPosName("NEventsPosMix_zBuf_");
327  TString hSibNegName("NEventsNegSib_zBuf_");
328  TString hMixNegName("NEventsNegMix_zBuf_");
329  while (zBin<99) {
330  TString fullSibName(hSibName.Data()); fullSibName+=zBin;
331  TString fullMixName(hMixName.Data()); fullMixName+=zBin;
332  TString fullSibPosName(hSibPosName.Data()); fullSibPosName+=zBin;
333  TString fullMixPosName(hMixPosName.Data()); fullMixPosName+=zBin;
334  TString fullSibNegName(hSibNegName.Data()); fullSibNegName+=zBin;
335  TString fullMixNegName(hMixNegName.Data()); fullMixNegName+=zBin;
336  inFile->cd();
337  inFile->GetObject(fullSibName.Data(),sib);
338  inFile->GetObject(fullMixName.Data(),mix);
339  TH1D *sibPos; inFile->GetObject(fullSibPosName.Data(),sibPos);
340  TH1D *mixPos; inFile->GetObject(fullMixPosName.Data(),mixPos);
341  TH1D *sibNeg; inFile->GetObject(fullSibNegName.Data(),sibNeg);
342  TH1D *mixNeg; inFile->GetObject(fullMixNegName.Data(),mixNeg);
343  if (sib && mix) {
344  outFile->cd();
345  sib->Write();
346  mix->Write();
347  sibPos->Write();
348  mixPos->Write();
349  sibNeg->Write();
350  mixNeg->Write();
351  } else {
352  break;
353  }
354  TString ptPName;
355  TString ptMName;
356  TString etaPName;
357  TString etaMName;
358  TH1 *tmp, *ptP, *ptM, *etaP, *etaM;
359  for (int iSpecies=0;iSpecies<2;iSpecies++) {
360  inFile->cd();
361  for (int iPar=0;iPar<nParentDist[iSpecies];iPar++) {
362  int jPar = parentDist[iPar][iSpecies];
363  ptPName = "meanPtP_parentBin"; ptPName += jPar; ptPName += "_zBuf_"; ptPName += zBin;
364  ptMName = "meanPtM_parentBin"; ptMName += jPar; ptMName += "_zBuf_"; ptMName += zBin;
365  etaPName = "etaP_parentBin"; etaPName += jPar; etaPName += "_zBuf_"; etaPName += zBin;
366  etaMName = "etaM_parentBin"; etaMName += jPar; etaMName += "_zBuf_"; etaMName += zBin;
367  if (0 == iPar) {
368  inFile->GetObject(ptPName.Data(),ptP);
369  inFile->GetObject(ptMName.Data(),ptM);
370  inFile->GetObject(etaPName.Data(),etaP);
371  inFile->GetObject(etaMName.Data(),etaM);
372  } else {
373  inFile->GetObject(ptPName.Data(),tmp); ptP->Add(tmp);
374  inFile->GetObject(ptMName.Data(),tmp); ptM->Add(tmp);
375  inFile->GetObject(etaPName.Data(),tmp); etaP->Add(tmp);
376  inFile->GetObject(etaMName.Data(),tmp); etaM->Add(tmp);
377  }
378  }
379  outFile->cd();
380  if (ptP) {
381  ptPName = "meanPtP"; ptPName += Species[iSpecies]; ptPName += "_zBuf_"; ptPName += zBin;
382  ptP->SetName(ptPName.Data());
383  ptP->Write();
384  }
385  if (ptM) {
386  ptMName = "meanPtM"; ptMName += Species[iSpecies]; ptMName += "_zBuf_"; ptMName += zBin;
387  ptM->SetName(ptMName.Data());
388  ptM->Write();
389  }
390  if (etaP) {
391  etaPName = "etaP"; etaPName += Species[iSpecies]; etaPName += "_zBuf_"; etaPName += zBin;
392  etaP->SetName(etaPName.Data());
393  etaP->Write();
394  }
395  if (etaM) {
396  etaMName = "etaM"; etaMName += Species[iSpecies]; etaMName += "_zBuf_"; etaMName += zBin;
397  etaM->SetName(etaMName.Data());
398  etaM->Write();
399  }
400  }
401  inFile->cd();
402  zBin++;
403  }
404  }
405  outFile->Close();
406 };
407 //--------------------------------------------------------------------------
408 void StEStructHAdd::addCuts(const char* outfile, const char* infile,
409  int* nlist, int ntot, int parentDist[][2], int* nParentDist, int symmXX) {
410 
411  TFile* tf=new TFile(infile);
412  if (tf) {
413  addCuts(outfile,tf,nlist,ntot,parentDist,nParentDist,symmXX);
414  }
415  return;
416 };
417 //------------------------------------------------------------------------
418 void StEStructHAdd::symmetrizeXX(TH2 *hist) {
419  for (int ix=1;ix<=hist->GetNbinsX();ix++) {
420  for (int iy=ix;iy<=hist->GetNbinsY();iy++) {
421  double sum = hist->GetBinContent(ix,iy) + hist->GetBinContent(iy,ix);
422  double err = sqrt(pow(hist->GetBinError(ix,iy),2) + pow(hist->GetBinError(iy,ix),2));
423  hist->SetBinContent(ix,iy,sum);
424  hist->SetBinContent(iy,ix,sum);
425  hist->SetBinError(ix,iy,err);
426  hist->SetBinError(iy,ix,err);
427  }
428  }
429 
430 };
431 //------------------------------------------------------------------------
432 void StEStructHAdd::old_addDensities(const char* outfile, TFile* inFile) {
433 
434  TFile* outFile = new TFile(outfile,"RECREATE");
435 
436  // If pair density histograms are in input file add z and cut bins,
437  // form LS and US, take ratio of sibling to mixed and write to output file.
438  inFile->cd();
439  TH2 *sibDens; inFile->GetObject("SibppTPCMidTZ_cutBin_0_zBuf_0",sibDens);
440  TH2 *mixDens; inFile->GetObject("MixppTPCMidTZ_cutBin_0_zBuf_0",mixDens);
441  if(sibDens && mixDens) {
442  int nCut = 0;
443  int nZBuf = 0;
444  TString hSibName("SibppTPCMidTZ_cutBin_");
445  while (nCut<99) {
446  TString full(hSibName.Data()); full += nCut; full += "_zBuf_0";
447  inFile->GetObject(full.Data(),sibDens);
448  if (sibDens) {
449  nCut++;
450  } else {
451  break;
452  }
453  }
454  while (nZBuf<99) {
455  TString full(hSibName.Data()); full += "0_zBuf_"; full += nZBuf;
456  inFile->GetObject(full.Data(),sibDens);
457  if (sibDens) {
458  nZBuf++;
459  } else {
460  break;
461  }
462  }
463 
464  TH2D *ls[6], *us[6], *tmp;
465  TH2D *cLS[6], *cUS[6];
466  inFile->GetObject("SibppTPCMidTZ_cutBin_0_zBuf_0",tmp);
467  for (int i=0;i<6;i++) {
468  ls[i] = (TH2D *) tmp->Clone(); ls[i]->Clear();
469  us[i] = (TH2D *) tmp->Clone(); us[i]->Clear();
470  cLS[i] = (TH2D *) tmp->Clone();
471  cUS[i] = (TH2D *) tmp->Clone();
472  }
473  TString hDensName;
474  const char *typePP[] = {"SibppTPCMidTZ_cutBin_", "SibppTPCMidTZC_cutBin_", "SibppTPCMidTZNC_cutBin_", "MixppTPCMidTZ_cutBin_", "MixppTPCMidTZC_cutBin_", "MixppTPCMidTZNC_cutBin_"};
475  const char *typeMM[] = {"SibmmTPCMidTZ_cutBin_", "SibmmTPCMidTZC_cutBin_", "SibmmTPCMidTZNC_cutBin_", "MixmmTPCMidTZ_cutBin_", "MixmmTPCMidTZC_cutBin_", "MixmmTPCMidTZNC_cutBin_"};
476  const char *typePM[] = {"SibpmTPCMidTZ_cutBin_", "SibpmTPCMidTZC_cutBin_", "SibpmTPCMidTZNC_cutBin_", "MixpmTPCMidTZ_cutBin_", "MixpmTPCMidTZC_cutBin_", "MixpmTPCMidTZNC_cutBin_"};
477  const char *typeMP[] = {"SibmpTPCMidTZ_cutBin_", "SibmpTPCMidTZC_cutBin_", "SibmpTPCMidTZNC_cutBin_", "MixmpTPCMidTZ_cutBin_", "MixmpTPCMidTZC_cutBin_", "MixmpTPCMidTZNC_cutBin_"};
478  const char *pairPos[] = {"TPCMidTZ_LS_", "TPCMidTZC_LS_", "TPCMidTZNC_LS_", "TPCMidTZ_US_", "TPCMidTZC_US_", "TPCMidTZNC_US_"};
479  for (int iCut=0;iCut<nCut;iCut++) {
480  inFile->cd();
481  for (int i=0;i<6;i++) {
482  cLS[i]->Clear();
483  cUS[i]->Clear();
484  }
485  hDensName = pairPos[0]; hDensName += iCut; cLS[0]->SetName(hDensName.Data());
486  hDensName = pairPos[1]; hDensName += iCut; cLS[1]->SetName(hDensName.Data());
487  hDensName = pairPos[2]; hDensName += iCut; cLS[2]->SetName(hDensName.Data());
488  hDensName = pairPos[3]; hDensName += iCut; cUS[0]->SetName(hDensName.Data());
489  hDensName = pairPos[4]; hDensName += iCut; cUS[1]->SetName(hDensName.Data());
490  hDensName = pairPos[5]; hDensName += iCut; cUS[2]->SetName(hDensName.Data());
491  for (int zBuf=0;zBuf<nZBuf;zBuf++) {
492  for (int it=0;it<6;it++) {
493  hDensName = typePP[it]; hDensName += iCut; hDensName += "_zBuf_"; hDensName += zBuf;
494  inFile->GetObject(hDensName.Data(),tmp);
495  cLS[it]->Add(tmp);
496  ls[it]->Add(tmp);
497  hDensName = typeMM[it]; hDensName += iCut; hDensName += "_zBuf_"; hDensName += zBuf;
498  inFile->GetObject(hDensName.Data(),tmp);
499  cLS[it]->Add(tmp);
500  ls[it]->Add(tmp);
501 
502  hDensName = typePM[it]; hDensName += iCut; hDensName += "_zBuf_"; hDensName += zBuf;
503  inFile->GetObject(hDensName.Data(),tmp);
504  cUS[it]->Add(tmp);
505  us[it]->Add(tmp);
506  hDensName = typeMP[it]; hDensName += iCut; hDensName += "_zBuf_"; hDensName += zBuf;
507  inFile->GetObject(hDensName.Data(),tmp);
508  cUS[it]->Add(tmp);
509  us[it]->Add(tmp);
510  }
511  }
512  cLS[0]->Divide(cLS[3]);
513  cLS[1]->Divide(cLS[4]);
514  cLS[2]->Divide(cLS[5]);
515  cUS[0]->Divide(cUS[3]);
516  cUS[1]->Divide(cUS[4]);
517  cUS[2]->Divide(cUS[5]);
518  outFile->cd();
519  cLS[0]->Write();
520  cLS[1]->Write();
521  cLS[2]->Write();
522  cUS[0]->Write();
523  cUS[1]->Write();
524  cUS[2]->Write();
525  }
526  ls[0]->Divide(ls[3]); ls[0]->SetName("TPCMidTZ_LS");
527  ls[1]->Divide(ls[4]); ls[1]->SetName("TPCMidTZC_LS");
528  ls[2]->Divide(ls[5]); ls[2]->SetName("TPCMidTZNC_LS");
529  us[0]->Divide(us[3]); us[0]->SetName("TPCMidTZ_US");
530  us[1]->Divide(us[4]); us[1]->SetName("TPCMidTZC_US");
531  us[2]->Divide(us[5]); us[2]->SetName("TPCMidTZNC_US");
532  outFile->cd();
533  ls[0]->Write();
534  ls[1]->Write();
535  ls[2]->Write();
536  us[0]->Write();
537  us[1]->Write();
538  us[2]->Write();
539  for (int i=0;i<6;i++) {
540  delete ls[i];
541  delete us[i];
542  delete cLS[i];
543  delete cUS[i];
544  }
545  }
546  outFile->Close();
547  delete outFile;
548 }
549 void StEStructHAdd::addDensities(const char* outfile, TFile* inFile) {
550 
551 
552  // If pair density histograms are in input file add z and cut bins,
553  // form LS and US, take ratio of sibling to mixed and write to output file.
554  TH2 *sibDens; inFile->GetObject("SibppTPCMidTZ_cutBin_0_zBuf_0",sibDens);
555  TH2 *mixDens; inFile->GetObject("MixppTPCMidTZ_cutBin_0_zBuf_0",mixDens);
556  if (!sibDens || !mixDens) {
557  return;
558  }
559 
560  // Find number of cut bins and number of zBuffers.
561  int nCut = 0;
562  int nZBuf = 0;
563  TString hSibName("SibppTPCMidTZ_cutBin_");
564  while (nCut<99) {
565  TString full(hSibName.Data()); full += nCut; full += "_zBuf_0";
566  inFile->GetObject(full.Data(),sibDens);
567  if (sibDens) {
568  nCut++;
569  } else {
570  break;
571  }
572  }
573  while (nZBuf<99) {
574  TString full(hSibName.Data()); full += "0_zBuf_"; full += nZBuf;
575  inFile->GetObject(full.Data(),sibDens);
576  if (sibDens) {
577  nZBuf++;
578  } else {
579  break;
580  }
581  }
582 
583 
584  TH1D *tmp1D;
585  TH2D *tmp2D;
586 
587  TH1D *sep1D[2][2][4][2];
588  TH1D *Sum_sep1D[2][2][4][2];
589  inFile->GetObject("SibppTPCAvgTSep_cutBin_0_zBuf_0",tmp1D);
590  for (int irel=0;irel<2;irel++) {
591  for (int ich=0;ich<2;ich++) {
592  for (int ipos=0;ipos<4;ipos++) {
593  for (int idir=0;idir<2;idir++) {
594  Sum_sep1D[irel][ich][ipos][idir] = (TH1D *) tmp1D->Clone();
595  }
596  }
597  }
598  }
599 
600  TH1D *dpt1D[2][2][2][2];
601  TH1D *Sum_dpt1D[2][2][2][2];
602  inFile->GetObject("SibppTPCMidTdptP_cutBin_0_zBuf_0",tmp1D);
603  for (int irel=0;irel<2;irel++) {
604  for (int ich=0;ich<2;ich++) {
605  for (int ity=0;ity<2;ity++) {
606  for (int idir=0;idir<2;idir++) {
607  Sum_dpt1D[irel][ich][ity][idir] = (TH1D *) tmp1D->Clone();
608  }
609  }
610  }
611  }
612 
613  TH2D *sep2D[2][2][4];
614  TH2D *Sum_sep2D[2][2][4];
615  inFile->GetObject("SibppTPCAvgTZ_cutBin_0_zBuf_0",tmp2D);
616  for (int irel=0;irel<2;irel++) {
617  for (int ich=0;ich<2;ich++) {
618  for (int ipos=0;ipos<4;ipos++) {
619  Sum_sep2D[irel][ich][ipos] = (TH2D *) tmp2D->Clone();
620  }
621  }
622  }
623 
624  TH2D *cross[2][2][2];
625  TH2D *Sum_cross[2][2][2];
626  inFile->GetObject("SibppTPCMidTZC_cutBin_0_zBuf_0",tmp2D);
627  for (int irel=0;irel<2;irel++) {
628  for (int ich=0;ich<2;ich++) {
629  for (int icr=0;icr<2;icr++) {
630  Sum_cross[irel][ich][icr] = (TH2D *) tmp2D->Clone();
631  }
632  }
633  }
634 
635  TH2D *dpt2D[2][2][3];
636  TH2D *Sum_dpt2D[2][2][3];
637  inFile->GetObject("SibppTPCEntTdpt_cutBin_0_zBuf_0",tmp2D);
638  for (int irel=0;irel<2;irel++) {
639  for (int ich=0;ich<2;ich++) {
640  for (int ipos=0;ipos<3;ipos++) {
641  Sum_dpt2D[irel][ich][ipos] = (TH2D *) tmp2D->Clone();
642  }
643  }
644  }
645 
646  TString hName;
647  const char *pos[] = {"Ent", "Mid", "Exit", "Avg"};
648  const char *dir[] = {"T", "Z"};
649  const char *ty[] = {"P", "N"};
650  const char *crs[] = {"C", "NC"};
651 
652  // Loop over cut bins.
653  // For each cut bin we sum over zBuffers, calculate and write Sib/Mix for LS and US
654  // We also sum over cut bins, calculating Sib/Mix at end.
655  const char *rel[] = {"Sib", "Mix"};
656  const char *ch[] = {"pp", "mm", "pm", "mp"};
657  const char *chType[] = {"LS", "US"};
658 
659  TFile* outFile = new TFile(outfile,"RECREATE");
660 
661  for (int iCut=0;iCut<nCut;iCut++) {
662  // Reset histograms for this cut. Will sum over zBuffers.
663  inFile->cd();
664  for (int irel=0;irel<2;irel++) {
665  for (int ich=0;ich<2;ich++) {
666  for (int ipos=0;ipos<4;ipos++) {
667  inFile->GetObject("SibppTPCAvgTSep_cutBin_0_zBuf_0",tmp1D);
668  for (int idir=0;idir<2;idir++) {
669  sep1D[irel][ich][ipos][idir] = (TH1D *) tmp1D->Clone();
670  sep1D[irel][ich][ipos][idir]->Clear();
671  }
672  }
673  inFile->GetObject("SibppTPCMidTdptP_cutBin_0_zBuf_0",tmp1D);
674  for (int ity=0;ity<2;ity++) {
675  for (int idir=0;idir<2;idir++) {
676  dpt1D[irel][ich][ity][idir] = (TH1D *) tmp1D->Clone();
677  dpt1D[irel][ich][ity][idir]->Clear();
678  }
679  }
680  inFile->GetObject("SibppTPCAvgTZ_cutBin_0_zBuf_0",tmp2D);
681  for (int ipos=0;ipos<4;ipos++) {
682  sep2D[irel][ich][ipos] = (TH2D *) tmp2D->Clone();
683  sep2D[irel][ich][ipos]->Clear();
684  }
685  inFile->GetObject("SibppTPCMidTZC_cutBin_0_zBuf_0",tmp2D);
686  for (int icr=0;icr<2;icr++) {
687  cross[irel][ich][icr] = (TH2D *) tmp2D->Clone();
688  cross[irel][ich][icr]->Clear();
689  }
690  inFile->GetObject("SibppTPCEntTdpt_cutBin_0_zBuf_0",tmp2D);
691  for (int ipos=0;ipos<3;ipos++) {
692  dpt2D[irel][ich][ipos] = (TH2D *) tmp2D->Clone();
693  dpt2D[irel][ich][ipos]->Clear();
694  }
695  }
696  }
697 
698  for (int zBuf=0;zBuf<nZBuf;zBuf++) {
699  // First do 1D T, Z separations.
700  for (int ipos=0;ipos<4;ipos++) {
701  for (int idir=0;idir<2;idir++) {
702  for (int irel=0;irel<2;irel++) {
703  for (int ich=0;ich<4;ich++) {
704  hName = rel[irel]; hName += ch[ich]; hName += "TPC"; hName += pos[ipos]; hName += dir[idir];
705  hName += "Sep_cutBin_"; hName += iCut; hName += "_zBuf_"; hName += zBuf;
706  inFile->GetObject(hName.Data(),tmp1D);
707  sep1D[irel][ich/2][ipos][idir]->Add(tmp1D);
708  Sum_sep1D[irel][ich/2][ipos][idir]->Add(tmp1D);
709  }
710  }
711  }
712  }
713  // Now 1D T*dpt, Z*dpt for P and N curvature
714  for (int ity=0;ity<2;ity++) {
715  for (int idir=0;idir<2;idir++) {
716  for (int irel=0;irel<2;irel++) {
717  for (int ich=0;ich<4;ich++) {
718  hName = rel[irel]; hName += ch[ich]; hName += "TPCMid"; hName += dir[idir]; hName += "dpt"; hName += ty[ity];
719  hName += "_cutBin_"; hName += iCut; hName += "_zBuf_"; hName += zBuf;
720  inFile->GetObject(hName.Data(),tmp1D);
721  dpt1D[irel][ich/2][ity][idir]->Add(tmp1D);
722  Sum_dpt1D[irel][ich/2][ity][idir]->Add(tmp1D);
723  }
724  }
725  }
726  }
727  // Now for 2D T-Z separations.
728  for (int ipos=0;ipos<4;ipos++) {
729  for (int irel=0;irel<2;irel++) {
730  for (int ich=0;ich<4;ich++) {
731  hName = rel[irel]; hName += ch[ich]; hName += "TPC"; hName += pos[ipos];
732  hName += "TZ_cutBin_"; hName += iCut; hName += "_zBuf_"; hName += zBuf;
733  inFile->GetObject(hName.Data(),tmp2D);
734  sep2D[irel][ich/2][ipos]->Add(tmp2D);
735  Sum_sep2D[irel][ich/2][ipos]->Add(tmp2D);
736  }
737  }
738  }
739  // Now 2D T, Z - dpt separations
740  for (int ipos=0;ipos<3;ipos++) {
741  for (int irel=0;irel<2;irel++) {
742  for (int ich=0;ich<4;ich++) {
743  hName = rel[irel]; hName += ch[ich]; hName += "TPC"; hName += pos[ipos];
744  hName += "Tdpt_cutBin_"; hName += iCut; hName += "_zBuf_"; hName += zBuf;
745  inFile->GetObject(hName.Data(),tmp2D);
746  dpt2D[irel][ich/2][ipos]->Add(tmp2D);
747  Sum_dpt2D[irel][ich/2][ipos]->Add(tmp2D);
748  }
749  }
750  }
751  // Finally, 2D T-Z separations dependent on crossing or not.
752  for (int icr=0;icr<2;icr++) {
753  for (int irel=0;irel<2;irel++) {
754  for (int ich=0;ich<4;ich++) {
755  hName = rel[irel]; hName += ch[ich]; hName += "TPCMidTZ"; hName += crs[icr];
756  hName += "_cutBin_"; hName += iCut; hName += "_zBuf_"; hName += zBuf;
757  inFile->GetObject(hName.Data(),tmp2D);
758  cross[irel][ich/2][icr]->Add(tmp2D);
759  Sum_cross[irel][ich/2][icr]->Add(tmp2D);
760  }
761  }
762  }
763  }
764 
765  // Calculate sibling/mixed (normalized) and write.
766  // Also want to set the name appropriately for when we read these from a file.
767  outFile->cd();
768 
769  for (int ipos=0;ipos<4;ipos++) {
770  for (int idir=0;idir<2;idir++) {
771  for (int ich=0;ich<2;ich++) {
772  sep1D[0][ich][ipos][idir]->Scale(sep1D[1][ich][ipos][idir]->Integral()/sep1D[0][ich][ipos][idir]->Integral());
773  sep1D[0][ich][ipos][idir]->Divide(sep1D[1][ich][ipos][idir]);
774  hName = chType[ich]; hName += pos[ipos]; hName += dir[idir]; hName += "Sep_cutBin_"; hName += iCut;
775  sep1D[0][ich][ipos][idir]->SetName(hName.Data());
776  sep1D[0][ich][ipos][idir]->SetTitle(hName.Data());
777  sep1D[0][ich][ipos][idir]->Write();
778  }
779  }
780  }
781  for (int ity=0;ity<2;ity++) {
782  for (int idir=0;idir<2;idir++) {
783  for (int ich=0;ich<2;ich++) {
784  dpt1D[0][ich][ity][idir]->Scale(dpt1D[1][ich][ity][idir]->Integral()/dpt1D[0][ich][ity][idir]->Integral());
785  dpt1D[0][ich][ity][idir]->Divide(dpt1D[1][ich][ity][idir]);
786  hName = chType[ich]; hName += dir[idir]; hName += ty[ity]; hName += "_cutBin_"; hName += iCut;
787  dpt1D[0][ich][ity][idir]->SetName(hName.Data());
788  dpt1D[0][ich][ity][idir]->SetTitle(hName.Data());
789  dpt1D[0][ich][ity][idir]->Write();
790  }
791  }
792  }
793  for (int ipos=0;ipos<4;ipos++) {
794  for (int ich=0;ich<2;ich++) {
795  sep2D[0][ich][ipos]->Scale(sep2D[1][ich][ipos]->Integral()/sep2D[0][ich][ipos]->Integral());
796  sep2D[0][ich][ipos]->Divide(sep2D[1][ich][ipos]);
797  hName = chType[ich]; hName += pos[ipos]; hName += "TZ_cutBin_"; hName += iCut;
798  sep2D[0][ich][ipos]->SetName(hName.Data());
799  sep2D[0][ich][ipos]->SetTitle(hName.Data());
800  sep2D[0][ich][ipos]->Write();
801  }
802  }
803  for (int ipos=0;ipos<3;ipos++) {
804  for (int ich=0;ich<2;ich++) {
805  dpt2D[0][ich][ipos]->Scale(dpt2D[1][ich][ipos]->Integral()/dpt2D[0][ich][ipos]->Integral());
806  dpt2D[0][ich][ipos]->Divide(dpt2D[1][ich][ipos]);
807  hName = chType[ich]; hName += pos[ipos]; hName += "Tdpt_cutBin_"; hName += iCut;
808  dpt2D[0][ich][ipos]->SetName(hName.Data());
809  dpt2D[0][ich][ipos]->SetTitle(hName.Data());
810  dpt2D[0][ich][ipos]->Write();
811  }
812  }
813  for (int icr=0;icr<2;icr++) {
814  for (int ich=0;ich<2;ich++) {
815  cross[0][ich][icr]->Scale(cross[1][ich][icr]->Integral()/cross[0][ich][icr]->Integral());
816  cross[0][ich][icr]->Divide(cross[1][ich][icr]);
817  hName = chType[ich]; hName += "MidTZ"; hName += crs[icr]; hName += "_cutBin_"; hName += iCut;
818  cross[0][ich][icr]->SetName(hName.Data());
819  cross[0][ich][icr]->SetTitle(hName.Data());
820  cross[0][ich][icr]->Write();
821  }
822  }
823 
824  for (int irel=0;irel<2;irel++) {
825  for (int ipos=0;ipos<4;ipos++) {
826  for (int idir=0;idir<2;idir++) {
827  for (int ich=0;ich<2;ich++) {
828  delete sep1D[irel][ich][ipos][idir];
829  }
830  }
831  }
832  for (int ity=0;ity<2;ity++) {
833  for (int idir=0;idir<2;idir++) {
834  for (int ich=0;ich<2;ich++) {
835  delete dpt1D[irel][ich][ity][idir];
836  }
837  }
838  }
839  for (int ipos=0;ipos<4;ipos++) {
840  for (int ich=0;ich<2;ich++) {
841  delete sep2D[irel][ich][ipos];
842  }
843  }
844  for (int ipos=0;ipos<3;ipos++) {
845  for (int ich=0;ich<2;ich++) {
846  delete dpt2D[irel][ich][ipos];
847  }
848  }
849  for (int icr=0;icr<2;icr++) {
850  for (int ich=0;ich<2;ich++) {
851  delete cross[irel][ich][icr];
852  }
853  }
854  }
855  }
856 
857  // Finished loop over cut bins.
858  // Now calculate and write ratios for summed cut bins.
859  outFile->cd();
860 
861  for (int ipos=0;ipos<4;ipos++) {
862  for (int idir=0;idir<2;idir++) {
863  for (int ich=0;ich<2;ich++) {
864  Sum_sep1D[0][ich][ipos][idir]->Scale(Sum_sep1D[1][ich][ipos][idir]->Integral()/Sum_sep1D[0][ich][ipos][idir]->Integral());
865  Sum_sep1D[0][ich][ipos][idir]->Divide(Sum_sep1D[1][ich][ipos][idir]);
866  hName = chType[ich]; hName += pos[ipos]; hName += dir[idir]; hName += "Sep";
867  Sum_sep1D[0][ich][ipos][idir]->SetName(hName.Data());
868  Sum_sep1D[0][ich][ipos][idir]->SetTitle(hName.Data());
869  Sum_sep1D[0][ich][ipos][idir]->Write();
870  }
871  }
872  }
873  for (int ity=0;ity<2;ity++) {
874  for (int idir=0;idir<2;idir++) {
875  for (int ich=0;ich<2;ich++) {
876  Sum_dpt1D[0][ich][ity][idir]->Scale(Sum_dpt1D[1][ich][ity][idir]->Integral()/Sum_dpt1D[0][ich][ity][idir]->Integral());
877  Sum_dpt1D[0][ich][ity][idir]->Divide(Sum_dpt1D[1][ich][ity][idir]);
878  hName = chType[ich]; hName += dir[idir]; hName += ty[ity];
879  Sum_dpt1D[0][ich][ity][idir]->SetName(hName.Data());
880  Sum_dpt1D[0][ich][ity][idir]->SetTitle(hName.Data());
881  Sum_dpt1D[0][ich][ity][idir]->Write();
882  }
883  }
884  }
885  for (int ipos=0;ipos<4;ipos++) {
886  for (int ich=0;ich<2;ich++) {
887  Sum_sep2D[0][ich][ipos]->Scale(Sum_sep2D[1][ich][ipos]->Integral()/Sum_sep2D[0][ich][ipos]->Integral());
888  Sum_sep2D[0][ich][ipos]->Divide(Sum_sep2D[1][ich][ipos]);
889  hName = chType[ich]; hName += pos[ipos]; hName += "TZ";
890  Sum_sep2D[0][ich][ipos]->SetName(hName.Data());
891  Sum_sep2D[0][ich][ipos]->SetTitle(hName.Data());
892  Sum_sep2D[0][ich][ipos]->Write();
893  }
894  }
895  for (int ipos=0;ipos<3;ipos++) {
896  for (int ich=0;ich<2;ich++) {
897  Sum_dpt2D[0][ich][ipos]->Scale(Sum_dpt2D[1][ich][ipos]->Integral()/Sum_dpt2D[0][ich][ipos]->Integral());
898  Sum_dpt2D[0][ich][ipos]->Divide(Sum_dpt2D[1][ich][ipos]);
899  hName = chType[ich]; hName += pos[ipos]; hName += "Tdpt";
900  Sum_dpt2D[0][ich][ipos]->SetName(hName.Data());
901  Sum_dpt2D[0][ich][ipos]->SetTitle(hName.Data());
902  Sum_dpt2D[0][ich][ipos]->Write();
903  }
904  }
905  for (int icr=0;icr<2;icr++) {
906  for (int ich=0;ich<2;ich++) {
907  Sum_cross[0][ich][icr]->Scale(Sum_cross[1][ich][icr]->Integral()/Sum_cross[0][ich][icr]->Integral());
908  Sum_cross[0][ich][icr]->Divide(Sum_cross[1][ich][icr]);
909  hName = chType[ich]; hName += "MidTZ"; hName += crs[icr];
910  Sum_cross[0][ich][icr]->SetName(hName.Data());
911  Sum_cross[0][ich][icr]->SetTitle(hName.Data());
912  Sum_cross[0][ich][icr]->Write();
913  }
914  }
915  // Maybe we need to delete histograms or garbage collection gets confused.
916  // (Maybe we just have to be careful to not confuse root into attempting to transfer
917  // histograms from input to output files.)
918 
919 
920  for (int irel=0;irel<2;irel++) {
921  for (int ipos=0;ipos<4;ipos++) {
922  for (int idir=0;idir<2;idir++) {
923  for (int ich=0;ich<2;ich++) {
924  delete Sum_sep1D[irel][ich][ipos][idir];
925  }
926  }
927  }
928  for (int ity=0;ity<2;ity++) {
929  for (int idir=0;idir<2;idir++) {
930  for (int ich=0;ich<2;ich++) {
931  delete Sum_dpt1D[irel][ich][ity][idir];
932  }
933  }
934  }
935  for (int ipos=0;ipos<4;ipos++) {
936  for (int ich=0;ich<2;ich++) {
937  delete Sum_sep2D[irel][ich][ipos];
938  }
939  }
940  for (int ipos=0;ipos<3;ipos++) {
941  for (int ich=0;ich<2;ich++) {
942  delete Sum_dpt2D[irel][ich][ipos];
943  }
944  }
945  for (int icr=0;icr<2;icr++) {
946  for (int ich=0;ich<2;ich++) {
947  delete Sum_cross[irel][ich][icr];
948  }
949  }
950  }
951 
952  outFile->Close();
953  delete outFile;
954 }
955 
956 //------------------------------------------------------------------------
957 void StEStructHAdd::combineUS(TFile * modFile) {
958 
959  // Add the -+ to the +- and clear the -+.
960  // This is for pid case when we are creating all. For non-identical particles we
961  // distinguish between +- and -+ but for identical particles we do not. In the
962  // case of 'all', which we make from summing all other cases, we don't
963  // want to distinguish +- and -+.
964 
965  const char* base[]={"Sib","Mix"};
966  const char* tpe[]={"pp","pm","mp","mm"};
967  const char* knd[]={"YtYt", "NYtYt", "PtPt",
968  "PhiPhi", "NPhiPhi", "PrPhiPhi", "PaPhiPhi", "PbPhiPhi",
969  "EtaEta", "PrEtaEta", "PaEtaEta", "PbEtaEta",
970  "DEtaDPhi", "NDEtaDPhi", "PrDEtaDPhi", "PaDEtaDPhi", "PbDEtaDPhi",
971  "SYtDYt", "NSYtDYt",
972  "SEtaDPhi", "NSEtaDPhi", "PrSEtaDPhi", "PaSEtaDPhi", "PbSEtaDPhi",
973  "Qinv", "NQinv"
974  "DEtaDPhiRot", "NDEtaDPhiRot", "PrDEtaDPhiRot", "PaDEtaDPhiRot", "PbDEtaDPhiRot"};
975 
976 // TFile* modFile = new TFile(fileName,"UPDATE");
977 
978  for (int i=0;i<2;i++) {
979  for (int k=0;k<28;k++) {
980 
981  TH1 *hpm=0, *hmp=0;
982  int zBin = 0;
983  while (zBin<99) {
984  if (k<31) {
985  // Histogram name is of form
986  // SibppDEtaDPhi_cutBin_n_zBuf_m
987  // where n is the usual cut bin number and m is the zbuffer index.
988  TString pmHist; pmHist+=base[i]; pmHist+=tpe[1]; pmHist+=knd[k];
989  pmHist+="_zBuf_"; pmHist+=zBin;
990  hpm=(TH1*)modFile->Get(pmHist.Data());
991  if (!hpm) {
992  // Have either processed all zbins or this histogram (at least for zBuf=0) isn't in the file.
993  goto lastZ;
994  }
995  TString mpHist; mpHist+=base[i]; mpHist+=tpe[2]; mpHist+=knd[k];
996  mpHist+="_zBuf_"; mpHist+=zBin;
997  hmp=(TH1*)modFile->Get(mpHist.Data());
998  hpm->Add(hmp);
999  hmp->Scale(0);
1000  modFile->cd();
1001  hpm->Write();
1002  hmp->Write();
1003  }
1004  zBin++;
1005  }
1006  lastZ:
1007  continue;
1008  }
1009  }
1010 // modFile->Close();
1011 };
1012 
1013 
1014 /***********************************************************************
1015  *
1016  * $Log: StEStructHAdd.cxx,v $
1017  * Revision 1.17 2012/11/16 21:27:23 prindle
1018  * 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.
1019  *
1020  * Revision 1.16 2010/06/23 22:33:41 prindle
1021  * In HAdd we distinguish between the parent distributions of the
1022  * two particles.
1023  * In Support I fixed a number of problems in the Pt correlation section.
1024  *
1025  * Revision 1.15 2010/03/02 21:48:30 prindle
1026  * Fix addDensities (for checking pair cuts)
1027  * Lots of small changes
1028  *
1029  * Revision 1.14 2009/11/09 21:32:59 prindle
1030  * Fix warnings about casting char * to a const char * by redeclaring as const char *.
1031  *
1032  * Revision 1.13 2009/05/08 00:21:42 prindle
1033  * In StEStructHadd remove support for old style of histogram names, do a better job calculating
1034  * errors (at least for number (\eta_\Delta,\phi_\Delta) histograms), double bins which
1035  * have an edge in the center (purely cosmetic when looking at intermediate histograms).
1036  * In StEStructSupport check for existance of histograms and return gracefully.
1037  * Code in buildChargeTypes and buildPtChargeTypes was essentially duplicate of code
1038  * in buildCommon and buildPtCommon so I refactored to reduce redundancy.
1039  *
1040  * Revision 1.12 2008/12/02 23:52:51 prindle
1041  * Get information about histogram XX being symmetrized from CutBin.
1042  * Changed TH1* to TH2D* in many places hoping to be able to plot DEtaDPhi
1043  * as colz (doesn't work yet).
1044  * Added direct calculation of \Delta\rho/\rho_{ref} (and similar) which is
1045  * needed for YtYt correlations.
1046  *
1047  * Revision 1.11 2008/05/01 23:46:39 prindle
1048  * Changed to use TH1D and TH2D (instead of TH1 and TH2) in some places so
1049  * we can use GetObject method to enforce type checking. Found I had missed
1050  * duplicating a \phi_\Delta row in one case. Also added a method to include
1051  * sum of pairdensity histograms in output file.
1052  *
1053  * Revision 1.10 2007/11/26 20:23:20 prindle
1054  * (Hadd to trick cvs to allow me to commit this modified file?)
1055  * Modifications to support saving of individual z-bins and the concept
1056  * of parent distributions important for pid analysis.
1057  *
1058  * Revision 1.9 2007/06/26 20:22:17 msd
1059  * Very minor bug fix
1060  *
1061  * Revision 1.8 2007/05/27 22:46:00 msd
1062  * Added buildChargeTypes mode 3 which takes rho_ref from track counts.
1063  * Added buildChargeTypeSumOfRatios.
1064  *
1065  * Revision 1.7 2007/01/26 17:20:57 msd
1066  * Updated HAdd for new binning scheme.
1067  * Improved Support::buildChargeTypes.
1068  *
1069  * Revision 1.6 2006/10/02 22:26:48 prindle
1070  * Hadd now symmetrizes histograms while adding them, so output is usable
1071  * in Support as before. Need to load library for Correlation so we know
1072  * how many bins there are.
1073  * Added alternative versions of methods to calculate Delta\sigma^2.
1074  * Important for pt correlations where we need proper normalization before
1075  * subtracting mixed reference.
1076  *
1077  * Revision 1.5 2006/04/06 01:09:46 prindle
1078  * Calculating pt for each cut bin caused changes in HAdd.
1079  * The splitting of +- into +- and -+ caused changes in Support.
1080  *
1081  * Revision 1.4 2005/09/29 17:42:14 msd
1082  * Added Xt to hadd
1083  *
1084  * Revision 1.3 2005/03/08 21:56:42 porter
1085  * fixed bug in StEStructHAdd.cxx and added diagnostic option in ptcorrelations to
1086  * view individual terms separately
1087  *
1088  * Revision 1.2 2005/03/03 01:33:04 porter
1089  * Added pt-correlations method to support and included
1090  * these histograms to the HAdd routine
1091  *
1092  * Revision 1.1 2004/07/01 00:37:13 porter
1093  * new code previously my StEStructHelper. Takes hists from correltation
1094  * pass and builds final ressults. Also the StEStructHAdd.h is a simple
1095  * replacemnt for my sumyt.C macro which could be expanded later as needed.
1096  *
1097  *
1098  *
1099  *********************************************************************/
1100 
1101 
1102