ATLAS Offline Software
LWHistRootUtils.icc
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 namespace LWHistRootUtilsInternals {
6 
7  template <class THX_LW, class THX>
8  inline THX* allocate1DRootHisto(THX_LW* lwhist,unsigned &arraysize) {
9  arraysize = lwhist->GetNbinsX()+2;
10  float * xbins = LWHistInt::getVarBins(lwhist);
11  return xbins
12  ? MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(), lwhist->GetNbinsX(), xbins)
13  : MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(), 1, lwhist->getXMin(), lwhist->getXMax());
14  }
15  template <class THX_LW, class THX>
16  inline THX* allocate2DRootHisto(THX_LW* lwhist,unsigned &arraysize) {
17  const unsigned nx(lwhist->GetNbinsX());
18  const unsigned ny(lwhist->GetNbinsY());
19  arraysize = (nx+2)*(ny+2);
20  const float * xbins_f = LWHistInt::getVarBinsX(lwhist);
21  const float * ybins_f = LWHistInt::getVarBinsY(lwhist);
22  double * xbins(0);
23  double * ybins(0);
24  if (xbins_f) {
25  xbins = LWPools::acquire<double>(nx+1);
26  for (unsigned i = 0;i<=nx;++i)
27  xbins[i]=xbins_f[i];
28  }
29  if (ybins_f) {
30  ybins = LWPools::acquire<double>(ny+1);
31  for (unsigned i = 0;i<=ny;++i)
32  ybins[i]=ybins_f[i];
33  }
34  THX * h(0);
35  if (ybins)
36  h = xbins
37  ? MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(),
38  lwhist->GetNbinsX(),xbins,
39  lwhist->GetNbinsY(),ybins)
40  : MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(),
41  lwhist->GetNbinsX(),lwhist->getXMin(),lwhist->getXMax(),
42  lwhist->GetNbinsY(),ybins);
43  else
44  h = xbins
45  ? MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(),
46  lwhist->GetNbinsX(),xbins,
47  lwhist->GetNbinsY(),lwhist->getYMin(),lwhist->getYMax())
48  : MP_NEW(THX)(lwhist->GetName(), lwhist->GetTitle(),
49  lwhist->GetNbinsX(),lwhist->getXMin(),lwhist->getXMax(),
50  lwhist->GetNbinsY(),lwhist->getYMin(),lwhist->getYMax());
51  assert(h);
52  if (xbins)
53  LWPools::release(xbins,nx+1);
54  if (ybins)
55  LWPools::release(ybins,ny+1);
56  return h;
57  //Fixme: don't use binlength of 1 for 2d??? Should speed up things...
58  }
59 
60  template <class THX_LW, class THX>
61  inline THX* allocateRootHisto(THX_LW* lwhist,unsigned &arraysize) { return allocate1DRootHisto<THX_LW,THX>(lwhist,arraysize); }
62  template <> inline TH2F* allocateRootHisto<TH2F_LW,TH2F>(TH2F_LW* lwhist,unsigned &arraysize) { return allocate2DRootHisto<TH2F_LW,TH2F>(lwhist,arraysize); }
63  template <> inline TH2D* allocateRootHisto<TH2D_LW,TH2D>(TH2D_LW* lwhist,unsigned &arraysize) { return allocate2DRootHisto<TH2D_LW,TH2D>(lwhist,arraysize); }
64  template <> inline TH2I* allocateRootHisto<TH2I_LW,TH2I>(TH2I_LW* lwhist,unsigned &arraysize) { return allocate2DRootHisto<TH2I_LW,TH2I>(lwhist,arraysize); }
65 
66  inline double * getArrayCopy(const TArrayD*a,unsigned &n)
67  {
68  n = a ? a->GetSize() : 0;
69  if (!n) return 0;
70  double * out = LWPools::acquire<double>(n);
71  for (unsigned i=0;i<n;++i)
72  out[i]=a->GetArray()[i];
73  return out;
74  }
75  inline void setHistParameters1D(TH1*hRoot,LWHist1D*lwhist) {
76  unsigned n;
77  double * xbins = getArrayCopy(hRoot->GetXaxis()->GetXbins(),n);
78  if (xbins) {
79  hRoot->SetBins(n-1, xbins);
80  LWPools::release(xbins,n);
81  } else {
82  hRoot->SetBins(lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax());
83  }
84  }
85  inline void setHistParameters2D(TH2*/*hRoot*/,LWHist2D*/*lwhist*/) {
86 //DISABLE// unsigned nx;
87 //DISABLE// double * xbins = getArrayCopy(hRoot->GetXaxis()->GetXbins(),nx);
88 //DISABLE// unsigned ny;
89 //DISABLE// double * ybins = getArrayCopy(hRoot->GetYaxis()->GetXbins(),ny);
90 //DISABLE// if (xbins) {
91 //DISABLE// if (ybins) {
92 //DISABLE// //Both X and Y have variable binning
93 //DISABLE// hRoot->SetBins(nx-1,xbins,
94 //DISABLE// ny-1,ybins);
95 //DISABLE// } else {
96 //DISABLE// //Just X (constructor exist, but not the corresponding setbins method)
97 //DISABLE// //assert(false);
98 //DISABLE// // hRoot->SetBins(nx-1,xbins,
99 //DISABLE// // lwhist->GetNbinsY(), lwhist->getYMin(), lwhist->getYMax());
100 //DISABLE// }
101 //DISABLE// } else if (ybins) {
102 //DISABLE// //Just Y (doesn't exist)
103 //DISABLE// // assert(false);
104 //DISABLE// // hRoot->SetBins(lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax(),
105 //DISABLE// // ny-1,ybins);
106 //DISABLE// } else {
107 //DISABLE// //None
108 //DISABLE// hRoot->SetBins(lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax(),
109 //DISABLE// lwhist->GetNbinsY(), lwhist->getYMin(), lwhist->getYMax());
110 //DISABLE// }
111 //DISABLE// if (xbins) LWPools::release(xbins,nx);
112 //DISABLE// if (ybins) LWPools::release(ybins,ny);
113  }
114  template <class THX_LW, class THX>
115  inline void setHistParameters(THX*roothist, THX_LW* lwhist) { setHistParameters1D(roothist,lwhist); }
116  template <> inline void setHistParameters<TH2F_LW,TH2F>(TH2F*roothist, TH2F_LW* lwhist) { setHistParameters2D(roothist,lwhist); }
117  template <> inline void setHistParameters<TH2D_LW,TH2D>(TH2D*roothist, TH2D_LW* lwhist) { setHistParameters2D(roothist,lwhist); }
118  template <> inline void setHistParameters<TH2I_LW,TH2I>(TH2I*roothist, TH2I_LW* lwhist) { setHistParameters2D(roothist,lwhist); }
119 
120  //Trick to avoid compile warnings:
121  template <class T> void multWithCast(T& destination, const double& fact) { destination = static_cast<T>(destination*fact); }
122 
123  template <class TArrayType>
124  void scaleArray(TArrayType&a, const double& fact) {
125  const unsigned n(a.fN);
126  if (n==0||!a.fArray)
127  return;
128  for (unsigned i=0;i<n;++i)
129  multWithCast(a.fArray[i],fact);//a.fArray[i] *= fact;
130  }
131 
132 }
133 
134 //____________________________________________________________________
135 template <class T, class TH_LW, class TH_root, class TFlexHist>
136 inline TH_root * LWHistRootUtils::createRootHisto(TH_LW* lwhist, TFlexHist * flexHist,bool& tookSumW2FromPools)
137 {
138  bool saveDefaultSumw2 = TH1::GetDefaultSumw2();
139  if (saveDefaultSumw2)
140  TH1::SetDefaultSumw2(false);
141 
142  unsigned arraysize;
143  TH_root * hRoot = LWHistRootUtilsInternals::allocateRootHisto<TH_LW,TH_root>(lwhist,arraysize);
144 
145  //Adopt() trick to allocate the array from our pools:
146  T * cont = LWPools::acquire<T>(arraysize);
147  tookSumW2FromPools = flexHist->holdsSeparateSumW2Info();
148  double * err = tookSumW2FromPools ? LWPools::acquire<double>(arraysize) : 0;
149 
150  flexHist->copyContents(cont,err);
151  hRoot->Adopt(arraysize,cont);
152  if (err) {
153  TArrayD& fSumw2 = LWHistRootUtils::getSumw2Array(hRoot);
154  assert(!fSumw2.fArray);
155  fSumw2.Adopt(arraysize,err);
156  }
157  LWHistRootUtilsInternals::setHistParameters(hRoot,lwhist);
158 
159  hRoot->SetEntries(flexHist->getEntries());
160  if (saveDefaultSumw2)
161  TH1::SetDefaultSumw2(true);
162 
163  return hRoot;
164 }
165 
166 //____________________________________________________________________
167 template <class THX>
168 void LWHistRootUtils::deleteRootHisto(THX*rootHist, bool& sumW2IsFromPools)
169 {
170  if (sumW2IsFromPools) {
171  TArrayD& fSumw2 = LWHistRootUtils::getSumw2Array(rootHist);
172  assert(fSumw2.fArray);
173  LWPools::release(fSumw2.fArray,rootHist->GetSize());
174  fSumw2.fArray = 0;
175  sumW2IsFromPools = false;
176  }
177  LWPools::release(rootHist->fArray,rootHist->GetSize());
178  rootHist->fArray=0;
179  MP_DELETE(rootHist);
180 }
181 
182 //____________________________________________________________________
183 inline TProfile * LWHistRootUtils::createRootProfileHisto(TProfile_LW* lwhist, Flex1DProfileHisto * flexHist)
184 {
185  const unsigned n(lwhist->GetNbinsX()+2);
186  double * entries = LWPools::acquire<double>(n);
187  double * contents = LWPools::acquire<double>(n);
188  double * errors = LWPools::acquire<double>(n);
189  flexHist->copyContents(entries,contents,errors);
190  float * xbins = LWHistInt::getVarBins(lwhist);
191  double ylow(flexHist->getProfParMin()), yup(flexHist->getProfParMax());
192  TProfile * hRoot(0);
193  const unsigned nbinsx(lwhist->GetNbinsX());
194  if (xbins) {
195  if (ylow==yup) {
196  hRoot=MP_NEW(TProfile)(lwhist->GetName(), lwhist->GetTitle(), nbinsx, xbins);
197  } else {
198  //Need double* array instead of float* array due to missing constructor:
199  double * xbins_d = LWPools::acquire<double>(nbinsx+1);
200  for (unsigned i = 0;i<=nbinsx;++i)
201  xbins_d[i]=xbins[i];
202  hRoot=MP_NEW(TProfile)(lwhist->GetName(), lwhist->GetTitle(), nbinsx, xbins_d,ylow,yup);
203  LWPools::release(xbins_d,nbinsx+1);
204  }
205  } else {
206  hRoot=MP_NEW(TProfile)(lwhist->GetName(), lwhist->GetTitle(), nbinsx, lwhist->getXMin(), lwhist->getXMax(),ylow,yup);
207  }
208  assert(hRoot);
209 
210  hRoot->Adopt(n,contents);
211  LWHistRootUtils::getSumw2Array(hRoot).Adopt(n,errors);
212  LWHistRootUtils::getBinEntriesArray(hRoot).Adopt(n,entries);
213  // hRoot->SetBins(lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax());
214  hRoot->SetEntries(lwhist->GetEntries());
215  hRoot->SetErrorOption(lwhist->GetErrorOption());
216  return hRoot;
217 }
218 
219 //____________________________________________________________________
220 inline TProfile2D * LWHistRootUtils::createRoot2DProfileHisto(TProfile2D_LW* lwhist, Flex2DProfileHisto * flexHist)
221 {
222  const unsigned nx(lwhist->GetNbinsX());
223  const unsigned ny(lwhist->GetNbinsY());
224  const unsigned n((nx+2)*(ny+2));
225  double * entries = LWPools::acquire<double>(n);
226  double * contents = LWPools::acquire<double>(n);
227  double * errors = LWPools::acquire<double>(n);
228  flexHist->copyContents(entries,contents,errors);
229 
230  const float * xbins_f = LWHistInt::getVarBinsX(lwhist);
231  const float * ybins_f = LWHistInt::getVarBinsY(lwhist);
232  double * xbins(0);
233  double * ybins(0);
234  if (xbins_f) {
235  xbins = LWPools::acquire<double>(nx+1);
236  for (unsigned i = 0;i<=nx;++i)
237  xbins[i]=xbins_f[i];
238  }
239  if (ybins_f) {
240  ybins = LWPools::acquire<double>(ny+1);
241  for (unsigned i = 0;i<=ny;++i)
242  ybins[i]=ybins_f[i];
243  }
244  TProfile2D * hRoot(0);
245 
246  if (ybins)
247  hRoot = xbins
248  ? MP_NEW(TProfile2D)( lwhist->GetName(), lwhist->GetTitle(),
249  lwhist->GetNbinsX(), xbins,
250  lwhist->GetNbinsY(), ybins )
251  : MP_NEW(TProfile2D)( lwhist->GetName(), lwhist->GetTitle(),
252  lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax(),
253  lwhist->GetNbinsY(), ybins );
254  else
255  hRoot = xbins
256  ? MP_NEW(TProfile2D)( lwhist->GetName(), lwhist->GetTitle(),
257  lwhist->GetNbinsX(), xbins,
258  lwhist->GetNbinsY(), lwhist->getYMin(), lwhist->getYMax() )
259  : MP_NEW(TProfile2D)( lwhist->GetName(), lwhist->GetTitle(),
260  lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax(),
261  lwhist->GetNbinsY(), lwhist->getYMin(), lwhist->getYMax(),
262  flexHist->getProfParMin(),flexHist->getProfParMax());
263 
264  assert(hRoot);
265  if (xbins)
266  LWPools::release(xbins,nx+1);
267  if (ybins)
268  LWPools::release(ybins,ny+1);
269  hRoot->Adopt(n,contents);
270  LWHistRootUtils::getSumw2Array(hRoot).Adopt(n,errors);
271  LWHistRootUtils::getBinEntriesArray(hRoot).Adopt(n,entries);
272  //hRoot->SetBins(lwhist->GetNbinsX(), lwhist->getXMin(), lwhist->getXMax(),
273  //lwhist->GetNbinsY(), lwhist->getYMin(), lwhist->getYMax());
274  hRoot->SetEntries(lwhist->GetEntries());
275  hRoot->SetErrorOption(lwhist->GetErrorOption());
276  return hRoot;
277 }
278 
279 //____________________________________________________________________
280 template <class TProfileX>
281 void LWHistRootUtils::deleteProfileHisto(TProfileX*rootHist)
282 {
283  const unsigned n(rootHist->GetSize());
284 
285  LWPools::release(getSumw2Array(rootHist).fArray,n);
286  getSumw2Array(rootHist).fArray = 0;
287 
288  LWPools::release(getBinEntriesArray(rootHist).fArray,n);
289  getBinEntriesArray(rootHist).fArray=0;
290 
291  LWPools::release(rootHist->fArray,n);
292  rootHist->fArray=0;
293 
294  MP_DELETE(rootHist);
295 }
296 
297 template <class THX>
298 void LWHistRootUtils::scaleContentsAndErrors( THX*h, const double& fact )
299 {
300  assert(h);
301  TArrayD& fSumw2 = LWHistRootUtils::getSumw2Array(h);
302  if (fSumw2.fN==0)
303  h->Sumw2();
304  LWHistRootUtilsInternals::scaleArray(fSumw2,fact*fact);
305  LWHistRootUtilsInternals::scaleArray(*h,fact);
306 }
Generated on Thu Jul 4 2024 21:21:04 for ATLAS Offline Software by   doxygen 1.8.18