Skip to content
This repository was archived by the owner on Jun 4, 2021. It is now read-only.

Commit ec05705

Browse files
kreczkokreczko
committed
Merge pull request #7 from kreczko/tunfold-root6
Fixed RooUnfoldTUnfold for ROOT6
2 parents 380cb61 + ebd9516 commit ec05705

File tree

3 files changed

+413
-5
lines changed

3 files changed

+413
-5
lines changed

‎GNUmakefile‎

Lines changed: 5 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -50,9 +50,9 @@ out := $(shell echo "$(ROOTSYS)/test/Makefile.arch not found - trying a basic Li
5050
ARCH = $(shell $(ROOTCONFIG) --arch)
5151
ROOTLIBS = $(shell $(ROOTCONFIG) --libs)
5252
CXXFLAGS = $(shell $(ROOTCONFIG) --cflags)
53-
CXX = g++
53+
CXX = clang++
5454
CXXFLAGS += -Wall -fPIC
55-
LD = g++
55+
LD = clang++
5656
LDFLAGS =
5757
SOFLAGS = -shared
5858
ObjSuf = o
@@ -118,10 +118,10 @@ HAVE_TUNFOLD = 1
118118
endif
119119
endif
120120

121-
#ifneq ($(HAVE_TUNFOLD),1)
121+
ifneq ($(HAVE_TUNFOLD),1)
122122
CPPFLAGS += -DNOTUNFOLD
123-
#EXCLUDE += RooUnfoldTUnfold.cxx RooUnfoldTUnfold.h
124-
#endif
123+
EXCLUDE += RooUnfoldTUnfold.cxx RooUnfoldTUnfold.h
124+
endif
125125

126126
# RooUnfoldDagostini is an interface to D'Agostini's implementation
127127
# of his algorithm: http://www.roma1.infn.it/~dagos/bayes_distr.txt .

‎src/RooUnfoldTUnfold.cxx‎

Lines changed: 273 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,273 @@
1+
//=====================================================================-*-C++-*-
2+
// File and Version Information:
3+
// $Id: RooUnfoldTUnfold.cxx 309 2011-10-10 20:40:14Z T.J.Adye $
4+
//
5+
// Description:
6+
// Unfolding class using TUnfold from ROOT to do the actual unfolding.
7+
//
8+
// Authors: Richard Claridge <[email protected]> & Tim Adye <[email protected]>
9+
//
10+
//==============================================================================
11+
12+
//____________________________________________________________
13+
/* BEGIN_HTML
14+
<p>Uses the unfolding method implemented in ROOT's <a href="http://root.cern.ch/root/html/TUnfold.html">TUnfold</a> class
15+
<p>Only included in ROOT versions 5.22 and higher
16+
<p>Only able to reconstruct 1 dimensional distributions
17+
<p>Can account for bin migration and smearing
18+
<p>Errors come as a full covariance matrix.
19+
<p>Will sometimes warn of "unlinked" bins. These are bins with 0 entries and do not effect the results of the unfolding
20+
<p>Regularisation parameter can be either optimised internally by plotting log10(chi2 squared) against log10(tau). The 'kink' in this curve is deemed the optimum tau value. This value can also be set manually (FixTau)
21+
<p>The latest version (TUnfold 15 in ROOT 2.27.04) will not handle plots with an additional underflow bin. As a result overflows must be turned off
22+
if v15 of TUnfold is used. ROOT versions 5.26 or below use v13 and so should be safe to use overflows.</ul>
23+
END_HTML */
24+
25+
/////////////////////////////////////////////////////////////
26+
27+
#include "RooUnfoldTUnfold.h"
28+
#include "RooUnfoldResponse.h"
29+
30+
#include <iostream>
31+
32+
#include "TH1.h"
33+
#include "TH2.h"
34+
#include "TVectorD.h"
35+
#include "TMatrixD.h"
36+
#include "TUnfold.h"
37+
#include "TGraph.h"
38+
39+
40+
using std::cout;
41+
using std::cerr;
42+
using std::endl;
43+
44+
ClassImp (RooUnfoldTUnfold);
45+
46+
RooUnfoldTUnfold::RooUnfoldTUnfold (const RooUnfoldTUnfold& rhs)
47+
: RooUnfold (rhs)
48+
{
49+
// Copy constructor.
50+
Init();
51+
CopyData (rhs);
52+
}
53+
54+
RooUnfoldTUnfold::RooUnfoldTUnfold (const RooUnfoldResponse* res, const TH1* meas, TUnfold::ERegMode reg,
55+
const char* name, const char* title)
56+
: RooUnfold (res, meas, name, title),_reg_method(reg)
57+
{
58+
// Constructor with response matrix object and measured unfolding input histogram.
59+
Init();
60+
}
61+
62+
void
63+
RooUnfoldTUnfold::Destroy()
64+
{
65+
delete _unf;
66+
}
67+
68+
RooUnfoldTUnfold*
69+
RooUnfoldTUnfold::Clone (const char* newname) const
70+
{
71+
//Clones object
72+
RooUnfoldTUnfold* unfold= new RooUnfoldTUnfold(*this);
73+
if (newname && strlen(newname)) unfold->SetName(newname);
74+
return unfold;
75+
}
76+
77+
void
78+
RooUnfoldTUnfold::Reset()
79+
{
80+
//Resets all values
81+
Destroy();
82+
Init();
83+
RooUnfold::Reset();
84+
}
85+
86+
87+
void
88+
RooUnfoldTUnfold::Assign (const RooUnfoldTUnfold& rhs)
89+
{
90+
RooUnfold::Assign (rhs);
91+
CopyData (rhs);
92+
}
93+
94+
void
95+
RooUnfoldTUnfold::CopyData (const RooUnfoldTUnfold& rhs)
96+
{
97+
tau_set=rhs.tau_set;
98+
_tau=rhs._tau;
99+
_reg_method=rhs._reg_method;
100+
}
101+
102+
103+
void
104+
RooUnfoldTUnfold::Init()
105+
{
106+
//Sets error matrix
107+
tau_set=false;
108+
_tau=0;
109+
_unf=0;
110+
GetSettings();
111+
}
112+
113+
TUnfold*
114+
RooUnfoldTUnfold::Impl()
115+
{
116+
return _unf;
117+
}
118+
119+
120+
void
121+
RooUnfoldTUnfold::Unfold()
122+
{
123+
/* Does the unfolding. Uses the optimal value of the unfolding parameter unless a value has already been set using FixTau*/
124+
125+
if (_nm<_nt) cerr << "Warning: fewer measured bins than truth bins. TUnfold may not work correctly." << endl;
126+
if (_haveCovMes) cerr << "Warning: TUnfold does not account for bin-bin correlations on measured input" << endl;
127+
128+
bool oldstat= TH1::AddDirectoryStatus();
129+
TH1::AddDirectory (kFALSE);
130+
TH1D* meas= HistNoOverflow (_meas, _overflow);
131+
TH2D* Hres=_res->HresponseNoOverflow();
132+
TH1::AddDirectory (oldstat);
133+
134+
// Add inefficiencies to measured overflow bin
135+
TVectorD tru= _res->Vtruth();
136+
for (int j= 1; j<=_nt; j++) {
137+
double ntru= 0.0;
138+
for (int i= 1; i<=_nm; i++) {
139+
ntru += Hres->GetBinContent(i,j);
140+
}
141+
Hres->SetBinContent (_nm+1, j, tru[j-1]-ntru);
142+
}
143+
144+
// Subtract fakes from measured distribution
145+
if (_res->FakeEntries()) {
146+
TVectorD fakes= _res->Vfakes();
147+
double fac= _res->Vmeasured().Sum();
148+
if (fac!=0.0) fac= Vmeasured().Sum() / fac;
149+
if (_verbose>=1) cout << "Subtract " << fac*fakes.Sum() << " fakes from measured distribution" << endl;
150+
for (int i= 1; i<=_nm; i++)
151+
meas->SetBinContent (i, meas->GetBinContent(i)-(fac*fakes[i-1]));
152+
}
153+
154+
int ndim= _meas->GetDimension();
155+
TUnfold::ERegMode reg= _reg_method;
156+
if (ndim == 2 || ndim == 3) reg= TUnfold::kRegModeNone; // set explicitly
157+
158+
_unf= new TUnfold(Hres,TUnfold::kHistMapOutputVert,reg);
159+
160+
if (ndim == 2) {
161+
int nx= _meas->GetNbinsX(), ny= _meas->GetNbinsY();
162+
_unf->RegularizeBins2D (0, 1, nx, nx, ny, _reg_method);
163+
} else if (ndim == 3) {
164+
int nx= _meas->GetNbinsX(), ny= _meas->GetNbinsY(), nz= _meas->GetNbinsZ(), nxy= nx*ny;
165+
for (int i= 0; i<nx; i++) {
166+
_unf->RegularizeBins2D ( i, nx, ny, nxy, nz, _reg_method);
167+
}
168+
for (int i= 0; i<ny; i++) {
169+
_unf->RegularizeBins2D ( nx*i, 1, nx, nxy, nz, _reg_method);
170+
}
171+
for (int i= 0; i<nz; i++) {
172+
_unf->RegularizeBins2D (nxy*i, 1, nx, nx, ny, _reg_method);
173+
}
174+
}
175+
176+
int nScan=30;
177+
// use automatic L-curve scan: start with taumin=taumax=0.0
178+
double tauMin=0.0;
179+
double tauMax=0.0;
180+
int iBest;
181+
TSpline *logTauX,*logTauY;
182+
TGraph *lCurve;
183+
// this method scans the parameter tau and finds the kink in the L curve
184+
// finally, the unfolding is done for the best choice of tau
185+
#if ROOT_VERSION_CODE >= ROOT_VERSION(5,23,0) /* TUnfold v6 (included in ROOT 5.22) didn't have setInput return value */
186+
int stat= _unf->SetInput(meas);
187+
if(stat>=10000) {
188+
cerr<<"Unfolding result may be wrong: " << stat/10000 << " unconstrained output bins\n";
189+
}
190+
#else
191+
_unf->SetInput(meas);
192+
#endif
193+
//_unf->SetConstraint(TUnfold::kEConstraintArea);
194+
if (!tau_set){
195+
iBest=_unf->ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
196+
_tau=_unf->GetTau(); // save value, even if we don't use it unless tau_set
197+
cout <<"Lcurve scan chose tau= "<<_tau<<endl;
198+
}
199+
else{
200+
_unf->DoUnfold(_tau);
201+
}
202+
TH1* reco= new TH1F();
203+
reco->SetNameTitle("_rec","reconstructed dist");
204+
_unf->GetOutput(reco);
205+
_rec.ResizeTo (_nt);
206+
for (int i=0;i<_nt;i++){
207+
_rec(i)=(reco->GetBinContent(i+1));
208+
}
209+
delete meas;
210+
delete Hres;
211+
delete reco;
212+
_unfolded= true;
213+
_haveCov= false;
214+
}
215+
216+
void
217+
RooUnfoldTUnfold::GetCov()
218+
{
219+
//Gets Covariance matrix
220+
if (!_unf) return;
221+
TH2D* ematrix= new TH2D();
222+
ematrix->SetNameTitle("ematrix","error matrix");
223+
_unf->GetEmatrix(ematrix);
224+
_cov.ResizeTo (_nt,_nt);
225+
for (int i= 0; i<_nt; i++) {
226+
for (int j= 0; j<_nt; j++) {
227+
_cov(i,j)= ematrix->GetBinContent(i+1,j+1);
228+
}
229+
}
230+
delete ematrix;
231+
_haveCov= true;
232+
}
233+
234+
235+
void
236+
RooUnfoldTUnfold::FixTau(double tau)
237+
{
238+
// Fix regularisation parameter to a specified value
239+
_tau=tau;
240+
tau_set=true;
241+
}
242+
243+
void
244+
RooUnfoldTUnfold::SetRegMethod(TUnfold::ERegMode regmethod)
245+
{
246+
/*
247+
Specifies the regularisation method:
248+
249+
regemthod setting regularisation
250+
=========================== ==============
251+
TUnfold::kRegModeNone none
252+
TUnfold::kRegModeSize minimize the size of (x-x0)
253+
TUnfold::kRegModeDerivative minimize the 1st derivative of (x-x0)
254+
TUnfold::kRegModeCurvature minimize the 2nd derivative of (x-x0)
255+
*/
256+
_reg_method=regmethod;
257+
}
258+
259+
void
260+
RooUnfoldTUnfold::OptimiseTau()
261+
{
262+
// Choose optimal regularisation parameter
263+
tau_set=false;
264+
}
265+
266+
void
267+
RooUnfoldTUnfold::GetSettings()
268+
{
269+
_minparm=0;
270+
_maxparm=1;
271+
_stepsizeparm=1e-2;
272+
_defaultparm=2;
273+
}

0 commit comments

Comments
 (0)