Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TSVDUnfold.cxx
Go to the documentation of this file.
1// Author: Kerstin Tackmann, Andreas Hoecker, Heiko Lacker
2
3/**********************************************************************************
4 * *
5 * Project: TSVDUnfold - data unfolding based on Singular Value Decomposition *
6 * Package: ROOT *
7 * Class : TSVDUnfold *
8 * *
9 * Description: *
10 * Single class implementation of SVD data unfolding based on: *
11 * A. Hoecker, V. Kartvelishvili, *
12 * "SVD approach to data unfolding" *
13 * NIM A372, 469 (1996) [hep-ph/9509307] *
14 * *
15 * Authors: *
16 * Kerstin Tackmann <Kerstin.Tackmann@cern.ch> - CERN, Switzerland *
17 * Andreas Hoecker <Andreas.Hoecker@cern.ch> - CERN, Switzerland *
18 * Heiko Lacker <lacker@physik.hu-berlin.de> - Humboldt U, Germany *
19 * *
20 * Copyright (c) 2010: *
21 * CERN, Switzerland *
22 * Humboldt University, Germany *
23 * *
24 **********************************************************************************/
25
26////////////////////////////////////////////////////////////////////////////////
27
28/** \class TSVDUnfold
29 \ingroup Hist
30 SVD Approach to Data Unfolding
31 <p>
32 Reference: <a href="http://arXiv.org/abs/hep-ph/9509307">Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307]</a>
33 <p>
34 TSVDUnfold implements the singular value decomposition based unfolding method (see reference). Currently, the unfolding of one-dimensional histograms is supported, with the same number of bins for the measured and the unfolded spectrum.
35 <p>
36 The unfolding procedure is based on singular value decomposition of the response matrix. The regularisation of the unfolding is implemented via a discrete minimum-curvature condition.
37 <p>
38 Monte Carlo inputs:
39 <ul>
40 <li><tt>xini</tt>: true underlying spectrum (TH1D, n bins)
41 <li><tt>bini</tt>: reconstructed spectrum (TH1D, n bins)
42 <li><tt>Adet</tt>: response matrix (TH2D, nxn bins)
43 </ul>
44 Consider the unfolding of a measured spectrum <tt>bdat</tt> with covariance matrix <tt>Bcov</tt> (if not passed explicitly, a diagonal covariance will be built given the errors of <tt>bdat</tt>). The corresponding spectrum in the Monte Carlo is given by <tt>bini</tt>, with the true underlying spectrum given by <tt>xini</tt>. The detector response is described by <tt>Adet</tt>, with <tt>Adet</tt> filled with events (not probabilities) with the true observable on the y-axis and the reconstructed observable on the x-axis.
45 <p>
46 The measured distribution can be unfolded for any combination of resolution, efficiency and acceptance effects, provided an appropriate definition of <tt>xini</tt> and <tt>Adet</tt>.<br><br>
47 <p>
48 The unfolding can be performed by
49 \code{.cpp}
50 TSVDUnfold *tsvdunf = new TSVDUnfold( bdat, Bcov, bini, xini, Adet );
51 TH1D* unfresult = tsvdunf->Unfold( kreg );
52 \endcode
53 where <tt>kreg</tt> determines the regularisation of the unfolding. In general, overregularisation (too small <tt>kreg</tt>) will bias the unfolded spectrum towards the Monte Carlo input, while underregularisation (too large <tt>kreg</tt>) will lead to large fluctuations in the unfolded spectrum. The optimal regularisation can be determined following guidelines in <a href="http://arXiv.org/abs/hep-ph/9509307">Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307]</a> using the distribution of the <tt>|d_i|</tt> that can be obtained by <tt>tsvdunf->GetD()</tt> and/or using pseudo-experiments.
54 <p>
55 Covariance matrices on the measured spectrum (for either the total uncertainties or individual sources of uncertainties) can be propagated to covariance matrices using the <tt>GetUnfoldCovMatrix</tt> method, which uses pseudo experiments for the propagation. In addition, <tt>GetAdetCovMatrix</tt> allows for the propagation of the statistical uncertainties on the response matrix using pseudo experiments. The covariance matrix corresponding to <tt>Bcov</tt> is also computed as described in <a href="http://arXiv.org/abs/hep-ph/9509307">Nucl. Instrum. Meth. A372, 469 (1996) [hep-ph/9509307]</a> and can be obtained from <tt>tsvdunf->GetXtau()</tt> and its (regularisation independent) inverse from <tt>tsvdunf->GetXinv()</tt>. The distribution of singular values can be retrieved using <tt>tsvdunf->GetSV()</tt>.
56 <p>
57 See also the tutorial for a toy example.
58*/
59
60#include <iostream>
61
62#include "TSVDUnfold.h"
63#include "TH1D.h"
64#include "TH2D.h"
65#include "TDecompSVD.h"
66#include "TRandom3.h"
67#include "TMath.h"
68
70
71////////////////////////////////////////////////////////////////////////////////
72/// Alternative constructor
73/// User provides data and MC test spectra, as well as detector response matrix, diagonal covariance matrix of measured spectrum built from the uncertainties on measured spectrum
74
75TSVDUnfold::TSVDUnfold( const TH1D *bdat, const TH1D *bini, const TH1D *xini, const TH2D *Adet )
76: TObject (),
77fNdim (0),
78fDdim (2),
79fNormalize (kFALSE),
80fKReg (-1),
81fDHist (nullptr),
82fSVHist (nullptr),
83fXtau (nullptr),
84fXinv (nullptr),
85fBdat (bdat),
86fBini (bini),
87fXini (xini),
88fAdet (Adet),
89fToyhisto (nullptr),
90fToymat (nullptr),
91fToyMode (kFALSE),
92fMatToyMode (kFALSE)
93{
94 if (bdat->GetNbinsX() != bini->GetNbinsX() ||
95 bdat->GetNbinsX() != xini->GetNbinsX() ||
96 bdat->GetNbinsX() != Adet->GetNbinsX() ||
97 bdat->GetNbinsX() != Adet->GetNbinsY()) {
98 TString msg = "All histograms must have equal dimension.\n";
99 msg += Form( " Found: dim(bdat)=%i\n", bdat->GetNbinsX() );
100 msg += Form( " Found: dim(bini)=%i\n", bini->GetNbinsX() );
101 msg += Form( " Found: dim(xini)=%i\n", xini->GetNbinsX() );
102 msg += Form( " Found: dim(Adet)=%i,%i\n", Adet->GetNbinsX(), Adet->GetNbinsY() );
103 msg += "Please start again!";
104
105 Fatal( "Init", msg, "%s" );
106 }
107
108 fBcov = (TH2D*)fAdet->Clone("bcov");
109
110 for(int i=1; i<=fBdat->GetNbinsX(); i++){
112 for(int j=1; j<=fBdat->GetNbinsX(); j++){
113 if(i==j) continue;
114 fBcov->SetBinContent(i,j,0.);
115 }
116 }
117 // Get the input histos
118 fNdim = bdat->GetNbinsX();
119 fDdim = 2; // This is the derivative used to compute the curvature matrix
120}
121
122
123////////////////////////////////////////////////////////////////////////////////
124/// Default constructor
125/// Initialisation of TSVDUnfold
126/// User provides data and MC test spectra, as well as detector response matrix and the covariance matrix of the measured distribution
127
128TSVDUnfold::TSVDUnfold( const TH1D *bdat, TH2D* Bcov, const TH1D *bini, const TH1D *xini, const TH2D *Adet )
129: TObject (),
130fNdim (0),
131fDdim (2),
132fNormalize (kFALSE),
133fKReg (-1),
134fDHist (nullptr),
135fSVHist (nullptr),
136fXtau (nullptr),
137fXinv (nullptr),
138fBdat (bdat),
139fBcov (Bcov),
140fBini (bini),
141fXini (xini),
142fAdet (Adet),
143fToyhisto (nullptr),
144fToymat (nullptr),
145fToyMode (kFALSE),
146fMatToyMode (kFALSE)
147{
148 if (bdat->GetNbinsX() != bini->GetNbinsX() ||
149 bdat->GetNbinsX() != xini->GetNbinsX() ||
150 bdat->GetNbinsX() != Bcov->GetNbinsX() ||
151 bdat->GetNbinsX() != Bcov->GetNbinsY() ||
152 bdat->GetNbinsX() != Adet->GetNbinsX() ||
153 bdat->GetNbinsX() != Adet->GetNbinsY()) {
154 TString msg = "All histograms must have equal dimension.\n";
155 msg += Form( " Found: dim(bdat)=%i\n", bdat->GetNbinsX() );
156 msg += Form( " Found: dim(Bcov)=%i,%i\n", Bcov->GetNbinsX(), Bcov->GetNbinsY() );
157 msg += Form( " Found: dim(bini)=%i\n", bini->GetNbinsX() );
158 msg += Form( " Found: dim(xini)=%i\n", xini->GetNbinsX() );
159 msg += Form( " Found: dim(Adet)=%i,%i\n", Adet->GetNbinsX(), Adet->GetNbinsY() );
160 msg += "Please start again!";
161
162 Fatal( "Init", msg, "%s" );
163 }
164
165 // Get the input histos
166 fNdim = bdat->GetNbinsX();
167 fDdim = 2; // This is the derivative used to compute the curvature matrix
168}
169
170////////////////////////////////////////////////////////////////////////////////
171/// Copy constructor
172
174: TObject ( other ),
175fNdim (other.fNdim),
176fDdim (other.fDdim),
177fNormalize (other.fNormalize),
178fKReg (other.fKReg),
179fDHist (other.fDHist),
180fSVHist (other.fSVHist),
181fXtau (other.fXtau),
182fXinv (other.fXinv),
183fBdat (other.fBdat),
184fBcov (other.fBcov),
185fBini (other.fBini),
186fXini (other.fXini),
187fAdet (other.fAdet),
188fToyhisto (other.fToyhisto),
189fToymat (other.fToymat),
190fToyMode (other.fToyMode),
191fMatToyMode (other.fMatToyMode)
192{
193}
194
195////////////////////////////////////////////////////////////////////////////////
196/// Destructor
197
199{
200 if(fToyhisto){
201 delete fToyhisto;
202 fToyhisto = nullptr;
203 }
204
205 if(fToymat){
206 delete fToymat;
207 fToymat = nullptr;
208 }
209
210 if(fDHist){
211 delete fDHist;
212 fDHist = nullptr;
213 }
214
215 if(fSVHist){
216 delete fSVHist;
217 fSVHist = nullptr;
218 }
219
220 if(fXtau){
221 delete fXtau;
222 fXtau = nullptr;
223 }
224
225 if(fXinv){
226 delete fXinv;
227 fXinv = nullptr;
228 }
229
230 if(fBcov){
231 delete fBcov;
232 fBcov = nullptr;
233 }
234}
235
236////////////////////////////////////////////////////////////////////////////////
237/// Perform the unfolding with regularisation parameter kreg
238
240{
241 fKReg = kreg;
242
243 // Make the histos
244 if (!fToyMode && !fMatToyMode) InitHistos( );
245
246 // Create vectors and matrices
247 TVectorD vb(fNdim), vbini(fNdim), vxini(fNdim), vberr(fNdim);
248 TMatrixD mB(fNdim, fNdim), mA(fNdim, fNdim), mCurv(fNdim, fNdim), mC(fNdim, fNdim);
249
250 Double_t eps = 1e-12;
251 Double_t sreg;
252
253 // Copy histograms entries into vector
254 if (fToyMode) { H2V( fToyhisto, vb ); H2Verr( fToyhisto, vberr ); }
255 else { H2V( fBdat, vb ); H2Verr( fBdat, vberr ); }
256
257 H2M( fBcov, mB);
258 H2V( fBini, vbini );
259 H2V( fXini, vxini );
260 if (fMatToyMode) H2M( fToymat, mA );
261 else H2M( fAdet, mA );
262
263 // Fill and invert the second derivative matrix
264 FillCurvatureMatrix( mCurv, mC );
265
266 // Inversion of mC by help of SVD
267 TDecompSVD CSVD(mC);
268 TMatrixD CUort = CSVD.GetU();
269 TMatrixD CVort = CSVD.GetV();
270 TVectorD CSV = CSVD.GetSig();
271
272 TMatrixD CSVM(fNdim, fNdim);
273 for (Int_t i=0; i<fNdim; i++) CSVM(i,i) = 1/CSV(i);
274
275 CUort.Transpose( CUort );
276 TMatrixD mCinv = (CVort*CSVM)*CUort;
277
278 //Rescale using the data covariance matrix
279 TDecompSVD BSVD( mB );
280 TMatrixD QT = BSVD.GetU();
281 QT.Transpose(QT);
282 TVectorD B2SV = BSVD.GetSig();
283 TVectorD BSV(B2SV);
284
285 for(int i=0; i<fNdim; i++){
286 BSV(i) = TMath::Sqrt(B2SV(i));
287 }
288 TMatrixD mAtmp(fNdim,fNdim);
289 TVectorD vbtmp(fNdim);
290 mAtmp *= 0;
291 vbtmp *= 0;
292 for(int i=0; i<fNdim; i++){
293 for(int j=0; j<fNdim; j++){
294 if(BSV(i)){
295 vbtmp(i) += QT(i,j)*vb(j)/BSV(i);
296 }
297 for(int m=0; m<fNdim; m++){
298 if(BSV(i)){
299 mAtmp(i,j) += QT(i,m)*mA(m,j)/BSV(i);
300 }
301 }
302 }
303 }
304 mA = mAtmp;
305 vb = vbtmp;
306
307 // Singular value decomposition and matrix operations
308 TDecompSVD ASVD( mA*mCinv );
309 TMatrixD Uort = ASVD.GetU();
310 TMatrixD Vort = ASVD.GetV();
311 TVectorD ASV = ASVD.GetSig();
312
313 if (!fToyMode && !fMatToyMode) {
314 V2H(ASV, *fSVHist);
315 }
316
317 TMatrixD Vreg = mCinv*Vort;
318
319 Uort.Transpose(Uort);
320 TVectorD vd = Uort*vb;
321
322 if (!fToyMode && !fMatToyMode) {
323 V2H(vd, *fDHist);
324 }
325
326 // Damping coefficient
327 Int_t k = GetKReg()-1;
328
329 TVectorD vx(fNdim); // Return variable
330
331 // Damping factors
332 TVectorD vdz(fNdim);
333 TMatrixD Z(fNdim, fNdim);
334 for (Int_t i=0; i<fNdim; i++) {
335 if (ASV(i)<ASV(0)*eps) sreg = ASV(0)*eps;
336 else sreg = ASV(i);
337 vdz(i) = sreg/(sreg*sreg + ASV(k)*ASV(k));
338 Z(i,i) = vdz(i)*vdz(i);
339 }
340 TVectorD vz = CompProd( vd, vdz );
341
342 TMatrixD VortT(Vort);
343 VortT.Transpose(VortT);
344 TMatrixD W = mCinv*Vort*Z*VortT*mCinv;
345
346 TMatrixD Xtau(fNdim, fNdim);
347 TMatrixD Xinv(fNdim, fNdim);
348 Xtau *= 0;
349 Xinv *= 0;
350 for (Int_t i=0; i<fNdim; i++) {
351 for (Int_t j=0; j<fNdim; j++) {
352 Xtau(i,j) = vxini(i) * vxini(j) * W(i,j);
353
354 double a=0;
355 for (Int_t m=0; m<fNdim; m++) {
356 a += mA(m,i)*mA(m,j);
357 }
358 if(vxini(i) && vxini(j))
359 Xinv(i,j) = a/vxini(i)/vxini(j);
360 }
361 }
362
363 // Compute the weights
364 TVectorD vw = Vreg*vz;
365
366 // Rescale by xini
367 vx = CompProd( vw, vxini );
368
369 if(fNormalize){ // Scale result to unit area
370 Double_t scale = vx.Sum();
371 if (scale > 0){
372 vx *= 1.0/scale;
373 Xtau *= 1./scale/scale;
374 Xinv *= scale*scale;
375 }
376 }
377
378 if (!fToyMode && !fMatToyMode) {
379 M2H(Xtau, *fXtau);
380 M2H(Xinv, *fXinv);
381 }
382
383 // Get Curvature and also chi2 in case of MC unfolding
384 if (!fToyMode && !fMatToyMode) {
385 Info( "Unfold", "Unfolding param: %i",k+1 );
386 Info( "Unfold", "Curvature of weight distribution: %f", GetCurvature( vw, mCurv ) );
387 }
388
389 TH1D* h = (TH1D*)fBdat->Clone("unfoldingresult");
390 for(int i=1; i<=fNdim; i++){
391 h->SetBinContent(i,0.);
392 h->SetBinError(i,0.);
393 }
394 V2H( vx, *h );
395
396 return h;
397}
398
399////////////////////////////////////////////////////////////////////////////////
400/// Determine for given input error matrix covariance matrix of unfolded
401/// spectrum from toy simulation given the passed covariance matrix on measured spectrum
402/// "cov" - covariance matrix on the measured spectrum, to be propagated
403/// "ntoys" - number of pseudo experiments used for the propagation
404/// "seed" - seed for pseudo experiments
405/// Note that this covariance matrix will contain effects of forced normalisation if spectrum is normalised to unit area.
406
408{
409 fToyMode = true;
410 TH1D* unfres = nullptr;
411 TH2D* unfcov = (TH2D*)fAdet->Clone("unfcovmat");
412 unfcov->SetTitle("Toy covariance matrix");
413 for(int i=1; i<=fNdim; i++)
414 for(int j=1; j<=fNdim; j++)
415 unfcov->SetBinContent(i,j,0.);
416
417 // Code for generation of toys (taken from RooResult and modified)
418 // Calculate the elements of the upper-triangular matrix L that
419 // gives Lt*L = C, where Lt is the transpose of L (the "square-root method")
420 TMatrixD L(fNdim,fNdim); L *= 0;
421
422 for (Int_t iPar= 0; iPar < fNdim; iPar++) {
423
424 // Calculate the diagonal term first
425 L(iPar,iPar) = cov->GetBinContent(iPar+1,iPar+1);
426 for (Int_t k=0; k<iPar; k++) L(iPar,iPar) -= TMath::Power( L(k,iPar), 2 );
427 if (L(iPar,iPar) > 0.0) L(iPar,iPar) = TMath::Sqrt(L(iPar,iPar));
428 else L(iPar,iPar) = 0.0;
429
430 // ...then the off-diagonal terms
431 for (Int_t jPar=iPar+1; jPar<fNdim; jPar++) {
432 L(iPar,jPar) = cov->GetBinContent(iPar+1,jPar+1);
433 for (Int_t k=0; k<iPar; k++) L(iPar,jPar) -= L(k,iPar)*L(k,jPar);
434 if (L(iPar,iPar)!=0.) L(iPar,jPar) /= L(iPar,iPar);
435 else L(iPar,jPar) = 0;
436 }
437 }
438
439 // Remember it
441 TRandom3 random(seed);
442
443 fToyhisto = (TH1D*)fBdat->Clone("toyhisto");
444 TH1D *toymean = (TH1D*)fBdat->Clone("toymean");
445 for (Int_t j=1; j<=fNdim; j++) toymean->SetBinContent(j,0.);
446
447 // Get the mean of the toys first
448 for (int i=1; i<=ntoys; i++) {
449
450 // create a vector of unit Gaussian variables
452 for (Int_t k= 0; k < fNdim; k++) g(k) = random.Gaus(0.,1.);
453
454 // Multiply this vector by Lt to introduce the appropriate correlations
455 g *= (*Lt);
456
457 // Add the mean value offsets and store the results
458 for (int j=1; j<=fNdim; j++) {
461 }
462
463 unfres = Unfold(GetKReg());
464
465 for (Int_t j=1; j<=fNdim; j++) {
466 toymean->SetBinContent(j, toymean->GetBinContent(j) + unfres->GetBinContent(j)/ntoys);
467 }
468 delete unfres;
469 unfres = nullptr;
470 }
471
472 // Reset the random seed
473 random.SetSeed(seed);
474
475 //Now the toys for the covariance matrix
476 for (int i=1; i<=ntoys; i++) {
477
478 // Create a vector of unit Gaussian variables
480 for (Int_t k= 0; k < fNdim; k++) g(k) = random.Gaus(0.,1.);
481
482 // Multiply this vector by Lt to introduce the appropriate correlations
483 g *= (*Lt);
484
485 // Add the mean value offsets and store the results
486 for (int j=1; j<=fNdim; j++) {
489 }
490 unfres = Unfold(GetKReg());
491
492 for (Int_t j=1; j<=fNdim; j++) {
493 for (Int_t k=1; k<=fNdim; k++) {
494 unfcov->SetBinContent(j,k,unfcov->GetBinContent(j,k) + ( (unfres->GetBinContent(j) - toymean->GetBinContent(j))* (unfres->GetBinContent(k) - toymean->GetBinContent(k))/(ntoys-1)) );
495 }
496 }
497 delete unfres;
498 unfres = nullptr;
499 }
500 delete Lt;
501 delete toymean;
503
504 return unfcov;
505}
506
507////////////////////////////////////////////////////////////////////////////////
508/// Determine covariance matrix of unfolded spectrum from finite statistics in
509/// response matrix using pseudo experiments
510/// "ntoys" - number of pseudo experiments used for the propagation
511/// "seed" - seed for pseudo experiments
512
514{
515 fMatToyMode = true;
516 TH1D* unfres = nullptr;
517 TH2D* unfcov = (TH2D*)fAdet->Clone("unfcovmat");
518 unfcov->SetTitle("Toy covariance matrix");
519 for(int i=1; i<=fNdim; i++)
520 for(int j=1; j<=fNdim; j++)
521 unfcov->SetBinContent(i,j,0.);
522
523 //Now the toys for the detector response matrix
524 TRandom3 random(seed);
525
526 fToymat = (TH2D*)fAdet->Clone("toymat");
527 TH1D *toymean = (TH1D*)fXini->Clone("toymean");
528 for (Int_t j=1; j<=fNdim; j++) toymean->SetBinContent(j,0.);
529
530 for (int i=1; i<=ntoys; i++) {
531 for (Int_t k=1; k<=fNdim; k++) {
532 for (Int_t m=1; m<=fNdim; m++) {
533 if (fAdet->GetBinContent(k,m)) {
535 }
536 }
537 }
538
539 unfres = Unfold(GetKReg());
540
541 for (Int_t j=1; j<=fNdim; j++) {
542 toymean->SetBinContent(j, toymean->GetBinContent(j) + unfres->GetBinContent(j)/ntoys);
543 }
544 delete unfres;
545 unfres = nullptr;
546 }
547
548 // Reset the random seed
549 random.SetSeed(seed);
550
551 for (int i=1; i<=ntoys; i++) {
552 for (Int_t k=1; k<=fNdim; k++) {
553 for (Int_t m=1; m<=fNdim; m++) {
554 if (fAdet->GetBinContent(k,m))
556 }
557 }
558
559 unfres = Unfold(GetKReg());
560
561 for (Int_t j=1; j<=fNdim; j++) {
562 for (Int_t k=1; k<=fNdim; k++) {
563 unfcov->SetBinContent(j,k,unfcov->GetBinContent(j,k) + ( (unfres->GetBinContent(j) - toymean->GetBinContent(j))*(unfres->GetBinContent(k) - toymean->GetBinContent(k))/(ntoys-1)) );
564 }
565 }
566 delete unfres;
567 unfres = nullptr;
568 }
569 delete toymean;
571
572 return unfcov;
573}
574
575////////////////////////////////////////////////////////////////////////////////
576/// Returns d vector (for choosing appropriate regularisation)
577
579{
580 for (int i=1; i<=fDHist->GetNbinsX(); i++) {
582 }
583 return fDHist;
584}
585
586////////////////////////////////////////////////////////////////////////////////
587/// Returns singular values vector
588
590{
591 return fSVHist;
592}
593
594////////////////////////////////////////////////////////////////////////////////
595/// Returns the computed regularized covariance matrix corresponding to total uncertainties on measured spectrum as passed in the constructor.
596/// Note that this covariance matrix will not contain the effects of forced normalization if spectrum is normalized to unit area.
597
599{
600 return fXtau;
601}
602
603////////////////////////////////////////////////////////////////////////////////
604/// Returns the computed inverse of the covariance matrix
605
607{
608 return fXinv;
609}
610
611////////////////////////////////////////////////////////////////////////////////
612/// Returns the covariance matrix
613
615{
616 return fBcov;
617}
618
619////////////////////////////////////////////////////////////////////////////////
620/// Fill 1D histogram into vector
621
622void TSVDUnfold::H2V( const TH1D* histo, TVectorD& vec )
623{
624 for (Int_t i=0; i<histo->GetNbinsX(); i++) vec(i) = histo->GetBinContent(i+1);
625}
626
627////////////////////////////////////////////////////////////////////////////////
628/// Fill 1D histogram errors into vector
629
630void TSVDUnfold::H2Verr( const TH1D* histo, TVectorD& vec )
631{
632 for (Int_t i=0; i<histo->GetNbinsX(); i++) vec(i) = histo->GetBinError(i+1);
633}
634
635////////////////////////////////////////////////////////////////////////////////
636/// Fill vector into 1D histogram
637
638void TSVDUnfold::V2H( const TVectorD& vec, TH1D& histo )
639{
640 for(Int_t i=0; i<vec.GetNrows(); i++) histo.SetBinContent(i+1, vec(i));
641}
642
643////////////////////////////////////////////////////////////////////////////////
644/// Fill 2D histogram into matrix
645
646void TSVDUnfold::H2M( const TH2D* histo, TMatrixD& mat )
647{
648 for (Int_t j=0; j<histo->GetNbinsX(); j++) {
649 for (Int_t i=0; i<histo->GetNbinsY(); i++) {
650 mat(i,j) = histo->GetBinContent(i+1,j+1);
651 }
652 }
653}
654
655////////////////////////////////////////////////////////////////////////////////
656/// Fill 2D histogram into matrix
657
658void TSVDUnfold::M2H( const TMatrixD& mat, TH2D& histo )
659{
660 for (Int_t j=0; j<mat.GetNcols(); j++) {
661 for (Int_t i=0; i<mat.GetNrows(); i++) {
662 histo.SetBinContent(i+1,j+1, mat(i,j));
663 }
664 }
665}
666
667////////////////////////////////////////////////////////////////////////////////
668/// Divide entries of two vectors
669
670TVectorD TSVDUnfold::VecDiv( const TVectorD& vec1, const TVectorD& vec2, Int_t zero )
671{
672 TVectorD quot(vec1.GetNrows());
673 for (Int_t i=0; i<vec1.GetNrows(); i++) {
674 if (vec2(i) != 0) quot(i) = vec1(i) / vec2(i);
675 else {
676 if (zero) quot(i) = 0;
677 else quot(i) = vec1(i);
678 }
679 }
680 return quot;
681}
682
683////////////////////////////////////////////////////////////////////////////////
684/// Divide matrix entries by vector
685
687{
688 TMatrixD quotmat(mat.GetNrows(), mat.GetNcols());
689 for (Int_t i=0; i<mat.GetNrows(); i++) {
690 for (Int_t j=0; j<mat.GetNcols(); j++) {
691 if (vec(i) != 0) quotmat(i,j) = mat(i,j) / vec(i);
692 else {
693 if (zero) quotmat(i,j) = 0;
694 else quotmat(i,j) = mat(i,j);
695 }
696 }
697 }
698 return quotmat;
699}
700
701////////////////////////////////////////////////////////////////////////////////
702/// Multiply entries of two vectors
703
705{
706 TVectorD res(vec1.GetNrows());
707 for (Int_t i=0; i<vec1.GetNrows(); i++) res(i) = vec1(i) * vec2(i);
708 return res;
709}
710
711////////////////////////////////////////////////////////////////////////////////
712/// Compute curvature of vector
713
715{
716 return vec*(curv*vec);
717}
718
719////////////////////////////////////////////////////////////////////////////////
720
722{
723 Double_t eps = 0.00001;
724
725 Int_t ndim = tCurv.GetNrows();
726
727 // Init
728 tCurv *= 0;
729 tC *= 0;
730
731 if (fDdim == 0) for (Int_t i=0; i<ndim; i++) tC(i,i) = 1;
732 else if (fDdim == 1) {
733 for (Int_t i=0; i<ndim; i++) {
734 if (i < ndim-1) tC(i,i+1) = 1.0;
735 tC(i,i) = 1.0;
736 }
737 }
738 else if (fDdim == 2) {
739 for (Int_t i=0; i<ndim; i++) {
740 if (i > 0) tC(i,i-1) = 1.0;
741 if (i < ndim-1) tC(i,i+1) = 1.0;
742 tC(i,i) = -2.0;
743 }
744 tC(0,0) = -1.0;
745 tC(ndim-1,ndim-1) = -1.0;
746 }
747 else if (fDdim == 3) {
748 for (Int_t i=1; i<ndim-2; i++) {
749 tC(i,i-1) = 1.0;
750 tC(i,i) = -3.0;
751 tC(i,i+1) = 3.0;
752 tC(i,i+2) = -1.0;
753 }
754 }
755 else if (fDdim==4) {
756 for (Int_t i=0; i<ndim; i++) {
757 if (i > 0) tC(i,i-1) = -4.0;
758 if (i < ndim-1) tC(i,i+1) = -4.0;
759 if (i > 1) tC(i,i-2) = 1.0;
760 if (i < ndim-2) tC(i,i+2) = 1.0;
761 tC(i,i) = 6.0;
762 }
763 tC(0,0) = 2.0;
764 tC(ndim-1,ndim-1) = 2.0;
765 tC(0,1) = -3.0;
766 tC(ndim-2,ndim-1) = -3.0;
767 tC(1,0) = -3.0;
768 tC(ndim-1,ndim-2) = -3.0;
769 tC(1,1) = 6.0;
770 tC(ndim-2,ndim-2) = 6.0;
771 }
772 else if (fDdim == 5) {
773 for (Int_t i=2; i < ndim-3; i++) {
774 tC(i,i-2) = 1.0;
775 tC(i,i-1) = -5.0;
776 tC(i,i) = 10.0;
777 tC(i,i+1) = -10.0;
778 tC(i,i+2) = 5.0;
779 tC(i,i+3) = -1.0;
780 }
781 }
782 else if (fDdim == 6) {
783 for (Int_t i = 3; i < ndim - 3; i++) {
784 tC(i,i-3) = 1.0;
785 tC(i,i-2) = -6.0;
786 tC(i,i-1) = 15.0;
787 tC(i,i) = -20.0;
788 tC(i,i+1) = 15.0;
789 tC(i,i+2) = -6.0;
790 tC(i,i+3) = 1.0;
791 }
792 }
793
794 // Add epsilon to avoid singularities
795 for (Int_t i=0; i<ndim; i++) tC(i,i) = tC(i,i) + eps;
796
797 //Get curvature matrix
798 for (Int_t i=0; i<ndim; i++) {
799 for (Int_t j=0; j<ndim; j++) {
800 for (Int_t k=0; k<ndim; k++) {
801 tCurv(i,j) = tCurv(i,j) + tC(k,i)*tC(k,j);
802 }
803 }
804 }
805}
806
807////////////////////////////////////////////////////////////////////////////////
808
810{
811 fDHist = new TH1D( "dd", "d vector after orthogonal transformation", fNdim, 0, fNdim );
812 fDHist->Sumw2();
813
814 fSVHist = new TH1D( "sv", "Singular values of AC^-1", fNdim, 0, fNdim );
815 fSVHist->Sumw2();
816
817 fXtau = (TH2D*)fAdet->Clone("Xtau");
818 fXtau->SetTitle("Regularized covariance matrix");
819 fXtau->Sumw2();
820
821 fXinv = (TH2D*)fAdet->Clone("Xinv");
822 fXinv->SetTitle("Inverse covariance matrix");
823 fXinv->Sumw2();
824}
825
826////////////////////////////////////////////////////////////////////////////////
827/// naive regularised inversion cuts off small elements
828
830{
831 // init reduced matrix
832 const UInt_t n = mat.GetNrows();
833 UInt_t nn = 0;
834
835 UInt_t *ipos = new UInt_t[n];
836 // UInt_t ipos[n];
837
838 // find max diagonal entries
839 Double_t ymax = 0;
840 for (UInt_t i=0; i<n; i++) if (TMath::Abs(mat[i][i]) > ymax) ymax = TMath::Abs(mat[i][i]);
841
842 for (UInt_t i=0; i<n; i++) {
843
844 // save position of accepted entries
845 if (TMath::Abs(mat[i][i])/ymax > eps) ipos[nn++] = i;
846 }
847
848 // effective matrix
849 TMatrixDSym matwork( nn );
850 for (UInt_t in=0; in<nn; in++) for (UInt_t jn=0; jn<nn; jn++) matwork(in,jn) = 0;
851
852 // fill non-zero effective working matrix
853 for (UInt_t in=0; in<nn; in++) {
854
855 matwork[in][in] = mat[ipos[in]][ipos[in]];
856 for (UInt_t jn=in+1; jn<nn; jn++) {
857 matwork[in][jn] = mat[ipos[in]][ipos[jn]];
858 matwork[jn][in] = matwork[in][jn];
859 }
860 }
861
862 // invert
863 matwork.Invert();
864
865 // reinitialise old matrix
866 for (UInt_t i=0; i<n; i++) for (UInt_t j=0; j<n; j++) mat[i][j] = 0;
867
868 // refill inverted matrix in old one
869 for (UInt_t in=0; in<nn; in++) {
870 mat[ipos[in]][ipos[in]] = matwork[in][in];
871 for (UInt_t jn=in+1; jn<nn; jn++) {
872 mat[ipos[in]][ipos[jn]] = matwork[in][jn];
873 mat[ipos[jn]][ipos[in]] = mat[ipos[in]][ipos[jn]];
874 }
875 }
876 delete [] ipos;
877}
878
879////////////////////////////////////////////////////////////////////////////////
880/// Helper routine to compute chi-squared between distributions using the computed inverse of the covariance matrix for the unfolded spectrum as given in paper.
881
882Double_t TSVDUnfold::ComputeChiSquared( const TH1D& truspec, const TH1D& unfspec)
883{
884 UInt_t n = truspec.GetNbinsX();
885
886 // compute chi2
887 Double_t chi2 = 0;
888 for (UInt_t i=0; i<n; i++) {
889 for (UInt_t j=0; j<n; j++) {
890 chi2 += ( (truspec.GetBinContent( i+1 )-unfspec.GetBinContent( i+1 )) *
891 (truspec.GetBinContent( j+1 )-unfspec.GetBinContent( j+1 )) * fXinv->GetBinContent(i+1,j+1) );
892 }
893 }
894
895 return chi2;
896}
897
#define g(i)
Definition RSha256.hxx:105
#define a(i)
Definition RSha256.hxx:99
#define h(i)
Definition RSha256.hxx:106
#define e(i)
Definition RSha256.hxx:103
constexpr Bool_t kFALSE
Definition RtypesCore.h:101
#define ClassImp(name)
Definition Rtypes.h:377
float ymax
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
char * Form(const char *fmt,...)
Formats a string in a circular formatting buffer.
Definition TString.cxx:2489
Single Value Decomposition class.
Definition TDecompSVD.h:24
const TVectorD & GetSig()
Definition TDecompSVD.h:59
const TMatrixD & GetV()
Definition TDecompSVD.h:57
const TMatrixD & GetU()
Definition TDecompSVD.h:55
1-D histogram with a double per channel (see TH1 documentation)
Definition TH1.h:669
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6686
virtual Int_t GetNbinsY() const
Definition TH1.h:298
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition TH1.cxx:9031
virtual Int_t GetNbinsX() const
Definition TH1.h:297
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition TH1.cxx:9174
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition TH1.cxx:9190
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5029
TObject * Clone(const char *newname="") const override
Make a complete copy of the underlying object.
Definition TH1.cxx:2752
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition TH1.cxx:8988
2-D histogram with a double per channel (see TH1 documentation)
Definition TH2.h:357
void SetBinContent(Int_t bin, Double_t content) override
Set bin content.
Definition TH2.cxx:2616
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition TH2.h:93
Int_t GetNrows() const
Int_t GetNcols() const
TMatrixTSym< Element > & Invert(Double_t *det=nullptr)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
Mother of all ROOT objects.
Definition TObject.h:41
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1015
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:961
Random number generator class based on M.
Definition TRandom3.h:27
void SetSeed(ULong_t seed=0) override
Set the random generator sequence if seed is 0 (default value) a TUUID is generated and used to fill ...
Definition TRandom3.cxx:206
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition TRandom.cxx:275
virtual ULong64_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition TRandom.cxx:404
SVD Approach to Data Unfolding.
Definition TSVDUnfold.h:46
TH1D * GetSV() const
Returns singular values vector.
TH2D * GetBCov() const
Returns the covariance matrix.
TH1D * fSVHist
! Distribution of singular values
Definition TSVDUnfold.h:133
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Bool_t fToyMode
! Internal switch for covariance matrix propagation
Definition TSVDUnfold.h:151
static void V2H(const TVectorD &vec, TH1D &histo)
Fill vector into 1D histogram.
static void H2M(const TH2D *histo, TMatrixD &mat)
Fill 2D histogram into matrix.
static void RegularisedSymMatInvert(TMatrixDSym &mat, Double_t eps=1e-3)
naive regularised inversion cuts off small elements
static TVectorD CompProd(const TVectorD &vec1, const TVectorD &vec2)
Multiply entries of two vectors.
static void M2H(const TMatrixD &mat, TH2D &histo)
Fill 2D histogram into matrix.
Bool_t fMatToyMode
! Internal switch for evaluation of statistical uncertainties from response matrix
Definition TSVDUnfold.h:152
TSVDUnfold(const TH1D *bdat, const TH1D *bini, const TH1D *xini, const TH2D *Adet)
Alternative constructor User provides data and MC test spectra, as well as detector response matrix,...
Bool_t fNormalize
! Normalize unfolded spectrum to 1
Definition TSVDUnfold.h:130
const TH1D * fBini
Reconstructed distribution (MC)
Definition TSVDUnfold.h:142
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Double_t ComputeChiSquared(const TH1D &truspec, const TH1D &unfspec)
Helper routine to compute chi-squared between distributions using the computed inverse of the covaria...
TH2D * fToymat
! Toy MC detector response matrix
Definition TSVDUnfold.h:150
static Double_t GetCurvature(const TVectorD &vec, const TMatrixD &curv)
Compute curvature of vector.
const TH1D * fBdat
Measured distribution (data)
Definition TSVDUnfold.h:140
~TSVDUnfold() override
Destructor.
static TMatrixD MatDivVec(const TMatrixD &mat, const TVectorD &vec, Int_t zero=0)
Divide matrix entries by vector.
TH2D * fXinv
! Computed inverse of covariance matrix
Definition TSVDUnfold.h:135
static TVectorD VecDiv(const TVectorD &vec1, const TVectorD &vec2, Int_t zero=0)
Divide entries of two vectors.
void FillCurvatureMatrix(TMatrixD &tCurv, TMatrixD &tC) const
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Int_t GetKReg() const
Definition TSVDUnfold.h:86
TH2D * fXtau
! Computed regularized covariance matrix
Definition TSVDUnfold.h:134
TH1D * fDHist
! Distribution of d (for checking regularization)
Definition TSVDUnfold.h:132
static void H2Verr(const TH1D *histo, TVectorD &vec)
Fill 1D histogram errors into vector.
Int_t fDdim
! Derivative for curvature matrix
Definition TSVDUnfold.h:129
TH2D * GetUnfoldCovMatrix(const TH2D *cov, Int_t ntoys, Int_t seed=1)
Determine for given input error matrix covariance matrix of unfolded spectrum from toy simulation giv...
TH2D * GetAdetCovMatrix(Int_t ntoys, Int_t seed=1)
Determine covariance matrix of unfolded spectrum from finite statistics in response matrix using pseu...
const TH2D * fAdet
Detector response matrix.
Definition TSVDUnfold.h:144
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
static void H2V(const TH1D *histo, TVectorD &vec)
Fill 1D histogram into vector.
Int_t fKReg
! Regularisation parameter
Definition TSVDUnfold.h:131
TH1D * fToyhisto
! Toy MC histogram
Definition TSVDUnfold.h:149
void InitHistos()
TH2D * fBcov
Covariance matrix of measured distribution (data)
Definition TSVDUnfold.h:141
const TH1D * fXini
Truth distribution (MC)
Definition TSVDUnfold.h:143
Int_t fNdim
! Truth and reconstructed dimensions
Definition TSVDUnfold.h:128
Basic string class.
Definition TString.h:139
Element Sum() const
Compute sum of elements.
Definition TVectorT.cxx:637
Int_t GetNrows() const
Definition TVectorT.h:73
const Int_t n
Definition legend1.C:16
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:662
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:721
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition TMathBase.h:123
TMarker m
Definition textangle.C:8