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 ) :
76fNdim (0),
77fDdim (2),
78fNormalize (kFALSE),
79fKReg (-1),
80fDHist (nullptr),
81fSVHist (nullptr),
82fXtau (nullptr),
83fXinv (nullptr),
84fBdat (bdat),
85fBini (bini),
86fXini (xini),
87fAdet (Adet),
88fToyhisto (nullptr),
89fToymat (nullptr),
90fToyMode (kFALSE),
91fMatToyMode (kFALSE)
92{
93 if (bdat->GetNbinsX() != bini->GetNbinsX() ||
94 bdat->GetNbinsX() != xini->GetNbinsX() ||
95 bdat->GetNbinsX() != Adet->GetNbinsX() ||
96 bdat->GetNbinsX() != Adet->GetNbinsY()) {
97 TString msg = "All histograms must have equal dimension.\n";
98 msg += Form( " Found: dim(bdat)=%i\n", bdat->GetNbinsX() );
99 msg += Form( " Found: dim(bini)=%i\n", bini->GetNbinsX() );
100 msg += Form( " Found: dim(xini)=%i\n", xini->GetNbinsX() );
101 msg += Form( " Found: dim(Adet)=%i,%i\n", Adet->GetNbinsX(), Adet->GetNbinsY() );
102 msg += "Please start again!";
103
104 Fatal( "Init", msg, "%s" );
105 }
106
107 fBcov = (TH2D*)fAdet->Clone("bcov");
108
109 for(int i=1; i<=fBdat->GetNbinsX(); i++){
111 for(int j=1; j<=fBdat->GetNbinsX(); j++){
112 if(i==j) continue;
113 fBcov->SetBinContent(i,j,0.);
114 }
115 }
116 // Get the input histos
117 fNdim = bdat->GetNbinsX();
118 fDdim = 2; // This is the derivative used to compute the curvature matrix
119}
120
121
122////////////////////////////////////////////////////////////////////////////////
123/// Default constructor
124/// Initialisation of TSVDUnfold
125/// User provides data and MC test spectra, as well as detector response matrix and the covariance matrix of the measured distribution
126
127TSVDUnfold::TSVDUnfold( const TH1D *bdat, TH2D* Bcov, const TH1D *bini, const TH1D *xini, const TH2D *Adet ) :
128fNdim (0),
129fDdim (2),
130fNormalize (kFALSE),
131fKReg (-1),
132fDHist (nullptr),
133fSVHist (nullptr),
134fXtau (nullptr),
135fXinv (nullptr),
136fBdat (bdat),
137fBcov (Bcov),
138fBini (bini),
139fXini (xini),
140fAdet (Adet),
141fToyhisto (nullptr),
142fToymat (nullptr),
143fToyMode (kFALSE),
144fMatToyMode (kFALSE)
145{
146 if (bdat->GetNbinsX() != bini->GetNbinsX() ||
147 bdat->GetNbinsX() != xini->GetNbinsX() ||
148 bdat->GetNbinsX() != Bcov->GetNbinsX() ||
149 bdat->GetNbinsX() != Bcov->GetNbinsY() ||
150 bdat->GetNbinsX() != Adet->GetNbinsX() ||
151 bdat->GetNbinsX() != Adet->GetNbinsY()) {
152 TString msg = "All histograms must have equal dimension.\n";
153 msg += Form( " Found: dim(bdat)=%i\n", bdat->GetNbinsX() );
154 msg += Form( " Found: dim(Bcov)=%i,%i\n", Bcov->GetNbinsX(), Bcov->GetNbinsY() );
155 msg += Form( " Found: dim(bini)=%i\n", bini->GetNbinsX() );
156 msg += Form( " Found: dim(xini)=%i\n", xini->GetNbinsX() );
157 msg += Form( " Found: dim(Adet)=%i,%i\n", Adet->GetNbinsX(), Adet->GetNbinsY() );
158 msg += "Please start again!";
159
160 Fatal( "Init", msg, "%s" );
161 }
162
163 // Get the input histos
164 fNdim = bdat->GetNbinsX();
165 fDdim = 2; // This is the derivative used to compute the curvature matrix
166}
167
168////////////////////////////////////////////////////////////////////////////////
169/// Copy constructor
170
172: TObject ( other ),
173fNdim (other.fNdim),
174fDdim (other.fDdim),
175fNormalize (other.fNormalize),
176fKReg (other.fKReg),
177fDHist (other.fDHist),
178fSVHist (other.fSVHist),
179fXtau (other.fXtau),
180fXinv (other.fXinv),
181fBdat (other.fBdat),
182fBcov (other.fBcov),
183fBini (other.fBini),
184fXini (other.fXini),
185fAdet (other.fAdet),
186fToyhisto (other.fToyhisto),
187fToymat (other.fToymat),
188fToyMode (other.fToyMode),
189fMatToyMode (other.fMatToyMode)
190{
191}
192
193////////////////////////////////////////////////////////////////////////////////
194/// Destructor
195
197{
198 if(fToyhisto){
199 delete fToyhisto;
200 fToyhisto = nullptr;
201 }
202
203 if(fToymat){
204 delete fToymat;
205 fToymat = nullptr;
206 }
207
208 if(fDHist){
209 delete fDHist;
210 fDHist = nullptr;
211 }
212
213 if(fSVHist){
214 delete fSVHist;
215 fSVHist = nullptr;
216 }
217
218 if(fXtau){
219 delete fXtau;
220 fXtau = nullptr;
221 }
222
223 if(fXinv){
224 delete fXinv;
225 fXinv = nullptr;
226 }
227
228 if(fBcov){
229 delete fBcov;
230 fBcov = nullptr;
231 }
232}
233
234////////////////////////////////////////////////////////////////////////////////
235/// Perform the unfolding with regularisation parameter kreg
236
238{
239 fKReg = kreg;
240
241 // Make the histos
242 if (!fToyMode && !fMatToyMode) InitHistos( );
243
244 // Create vectors and matrices
245 TVectorD vb(fNdim), vbini(fNdim), vxini(fNdim), vberr(fNdim);
246 TMatrixD mB(fNdim, fNdim), mA(fNdim, fNdim), mCurv(fNdim, fNdim), mC(fNdim, fNdim);
247
248 Double_t eps = 1e-12;
249 Double_t sreg;
250
251 // Copy histograms entries into vector
252 if (fToyMode) { H2V( fToyhisto, vb ); H2Verr( fToyhisto, vberr ); }
253 else { H2V( fBdat, vb ); H2Verr( fBdat, vberr ); }
254
255 H2M( fBcov, mB);
256 H2V( fBini, vbini );
257 H2V( fXini, vxini );
258 if (fMatToyMode) H2M( fToymat, mA );
259 else H2M( fAdet, mA );
260
261 // Fill and invert the second derivative matrix
262 FillCurvatureMatrix( mCurv, mC );
263
264 // Inversion of mC by help of SVD
265 TDecompSVD CSVD(mC);
266 TMatrixD CUort = CSVD.GetU();
267 TMatrixD CVort = CSVD.GetV();
268 TVectorD CSV = CSVD.GetSig();
269
270 TMatrixD CSVM(fNdim, fNdim);
271 for (Int_t i=0; i<fNdim; i++) CSVM(i,i) = 1/CSV(i);
272
273 CUort.Transpose( CUort );
274 TMatrixD mCinv = (CVort*CSVM)*CUort;
275
276 //Rescale using the data covariance matrix
277 TDecompSVD BSVD( mB );
278 TMatrixD QT = BSVD.GetU();
279 QT.Transpose(QT);
280 TVectorD B2SV = BSVD.GetSig();
281 TVectorD BSV(B2SV);
282
283 for(int i=0; i<fNdim; i++){
284 BSV(i) = TMath::Sqrt(B2SV(i));
285 }
286 TMatrixD mAtmp(fNdim,fNdim);
287 TVectorD vbtmp(fNdim);
288 mAtmp *= 0;
289 vbtmp *= 0;
290 for(int i=0; i<fNdim; i++){
291 for(int j=0; j<fNdim; j++){
292 if(BSV(i)){
293 vbtmp(i) += QT(i,j)*vb(j)/BSV(i);
294 }
295 for(int m=0; m<fNdim; m++){
296 if(BSV(i)){
297 mAtmp(i,j) += QT(i,m)*mA(m,j)/BSV(i);
298 }
299 }
300 }
301 }
302 mA = mAtmp;
303 vb = vbtmp;
304
305 // Singular value decomposition and matrix operations
306 TDecompSVD ASVD( mA*mCinv );
307 TMatrixD Uort = ASVD.GetU();
308 TMatrixD Vort = ASVD.GetV();
309 TVectorD ASV = ASVD.GetSig();
310
311 if (!fToyMode && !fMatToyMode) {
312 V2H(ASV, *fSVHist);
313 }
314
315 TMatrixD Vreg = mCinv*Vort;
316
317 Uort.Transpose(Uort);
318 TVectorD vd = Uort*vb;
319
320 if (!fToyMode && !fMatToyMode) {
321 V2H(vd, *fDHist);
322 }
323
324 // Damping coefficient
325 Int_t k = GetKReg()-1;
326
327 TVectorD vx(fNdim); // Return variable
328
329 // Damping factors
330 TVectorD vdz(fNdim);
331 TMatrixD Z(fNdim, fNdim);
332 for (Int_t i=0; i<fNdim; i++) {
333 if (ASV(i)<ASV(0)*eps) sreg = ASV(0)*eps;
334 else sreg = ASV(i);
335 vdz(i) = sreg/(sreg*sreg + ASV(k)*ASV(k));
336 Z(i,i) = vdz(i)*vdz(i);
337 }
338 TVectorD vz = CompProd( vd, vdz );
339
340 TMatrixD VortT(Vort);
341 VortT.Transpose(VortT);
342 TMatrixD W = mCinv*Vort*Z*VortT*mCinv;
343
344 TMatrixD Xtau(fNdim, fNdim);
345 TMatrixD Xinv(fNdim, fNdim);
346 Xtau *= 0;
347 Xinv *= 0;
348 for (Int_t i=0; i<fNdim; i++) {
349 for (Int_t j=0; j<fNdim; j++) {
350 Xtau(i,j) = vxini(i) * vxini(j) * W(i,j);
351
352 double a=0;
353 for (Int_t m=0; m<fNdim; m++) {
354 a += mA(m,i)*mA(m,j);
355 }
356 if(vxini(i) && vxini(j))
357 Xinv(i,j) = a/vxini(i)/vxini(j);
358 }
359 }
360
361 // Compute the weights
362 TVectorD vw = Vreg*vz;
363
364 // Rescale by xini
365 vx = CompProd( vw, vxini );
366
367 if(fNormalize){ // Scale result to unit area
368 Double_t scale = vx.Sum();
369 if (scale > 0){
370 vx *= 1.0/scale;
371 Xtau *= 1./scale/scale;
372 Xinv *= scale*scale;
373 }
374 }
375
376 if (!fToyMode && !fMatToyMode) {
377 M2H(Xtau, *fXtau);
378 M2H(Xinv, *fXinv);
379 }
380
381 // Get Curvature and also chi2 in case of MC unfolding
382 if (!fToyMode && !fMatToyMode) {
383 Info( "Unfold", "Unfolding param: %i",k+1 );
384 Info( "Unfold", "Curvature of weight distribution: %f", GetCurvature( vw, mCurv ) );
385 }
386
387 TH1D* h = (TH1D*)fBdat->Clone("unfoldingresult");
388 for(int i=1; i<=fNdim; i++){
389 h->SetBinContent(i,0.);
390 h->SetBinError(i,0.);
391 }
392 V2H( vx, *h );
393
394 return h;
395}
396
397////////////////////////////////////////////////////////////////////////////////
398/// Determine for given input error matrix covariance matrix of unfolded
399/// spectrum from toy simulation given the passed covariance matrix on measured spectrum
400/// "cov" - covariance matrix on the measured spectrum, to be propagated
401/// "ntoys" - number of pseudo experiments used for the propagation
402/// "seed" - seed for pseudo experiments
403/// Note that this covariance matrix will contain effects of forced normalisation if spectrum is normalised to unit area.
404
406{
407 fToyMode = true;
408 TH1D* unfres = nullptr;
409 TH2D* unfcov = (TH2D*)fAdet->Clone("unfcovmat");
410 unfcov->SetTitle("Toy covariance matrix");
411 for(int i=1; i<=fNdim; i++)
412 for(int j=1; j<=fNdim; j++)
413 unfcov->SetBinContent(i,j,0.);
414
415 // Code for generation of toys (taken from RooResult and modified)
416 // Calculate the elements of the upper-triangular matrix L that
417 // gives Lt*L = C, where Lt is the transpose of L (the "square-root method")
418 TMatrixD L(fNdim,fNdim); L *= 0;
419
420 for (Int_t iPar= 0; iPar < fNdim; iPar++) {
421
422 // Calculate the diagonal term first
423 L(iPar,iPar) = cov->GetBinContent(iPar+1,iPar+1);
424 for (Int_t k=0; k<iPar; k++) L(iPar,iPar) -= TMath::Power( L(k,iPar), 2 );
425 if (L(iPar,iPar) > 0.0) L(iPar,iPar) = TMath::Sqrt(L(iPar,iPar));
426 else L(iPar,iPar) = 0.0;
427
428 // ...then the off-diagonal terms
429 for (Int_t jPar=iPar+1; jPar<fNdim; jPar++) {
430 L(iPar,jPar) = cov->GetBinContent(iPar+1,jPar+1);
431 for (Int_t k=0; k<iPar; k++) L(iPar,jPar) -= L(k,iPar)*L(k,jPar);
432 if (L(iPar,iPar)!=0.) L(iPar,jPar) /= L(iPar,iPar);
433 else L(iPar,jPar) = 0;
434 }
435 }
436
437 // Remember it
439 TRandom3 random(seed);
440
441 fToyhisto = (TH1D*)fBdat->Clone("toyhisto");
442 TH1D *toymean = (TH1D*)fBdat->Clone("toymean");
443 for (Int_t j=1; j<=fNdim; j++) toymean->SetBinContent(j,0.);
444
445 // Get the mean of the toys first
446 for (int i=1; i<=ntoys; i++) {
447
448 // create a vector of unit Gaussian variables
450 for (Int_t k= 0; k < fNdim; k++) g(k) = random.Gaus(0.,1.);
451
452 // Multiply this vector by Lt to introduce the appropriate correlations
453 g *= (*Lt);
454
455 // Add the mean value offsets and store the results
456 for (int j=1; j<=fNdim; j++) {
459 }
460
461 unfres = Unfold(GetKReg());
462
463 for (Int_t j=1; j<=fNdim; j++) {
464 toymean->SetBinContent(j, toymean->GetBinContent(j) + unfres->GetBinContent(j)/ntoys);
465 }
466 delete unfres;
467 unfres = nullptr;
468 }
469
470 // Reset the random seed
471 random.SetSeed(seed);
472
473 //Now the toys for the covariance matrix
474 for (int i=1; i<=ntoys; i++) {
475
476 // Create a vector of unit Gaussian variables
478 for (Int_t k= 0; k < fNdim; k++) g(k) = random.Gaus(0.,1.);
479
480 // Multiply this vector by Lt to introduce the appropriate correlations
481 g *= (*Lt);
482
483 // Add the mean value offsets and store the results
484 for (int j=1; j<=fNdim; j++) {
487 }
488 unfres = Unfold(GetKReg());
489
490 for (Int_t j=1; j<=fNdim; j++) {
491 for (Int_t k=1; k<=fNdim; k++) {
492 unfcov->SetBinContent(j,k,unfcov->GetBinContent(j,k) + ( (unfres->GetBinContent(j) - toymean->GetBinContent(j))* (unfres->GetBinContent(k) - toymean->GetBinContent(k))/(ntoys-1)) );
493 }
494 }
495 delete unfres;
496 unfres = nullptr;
497 }
498 delete Lt;
499 delete toymean;
501
502 return unfcov;
503}
504
505////////////////////////////////////////////////////////////////////////////////
506/// Determine covariance matrix of unfolded spectrum from finite statistics in
507/// response matrix using pseudo experiments
508/// "ntoys" - number of pseudo experiments used for the propagation
509/// "seed" - seed for pseudo experiments
510
512{
513 fMatToyMode = true;
514 TH1D* unfres = nullptr;
515 TH2D* unfcov = (TH2D*)fAdet->Clone("unfcovmat");
516 unfcov->SetTitle("Toy covariance matrix");
517 for(int i=1; i<=fNdim; i++)
518 for(int j=1; j<=fNdim; j++)
519 unfcov->SetBinContent(i,j,0.);
520
521 //Now the toys for the detector response matrix
522 TRandom3 random(seed);
523
524 fToymat = (TH2D*)fAdet->Clone("toymat");
525 TH1D *toymean = (TH1D*)fXini->Clone("toymean");
526 for (Int_t j=1; j<=fNdim; j++) toymean->SetBinContent(j,0.);
527
528 for (int i=1; i<=ntoys; i++) {
529 for (Int_t k=1; k<=fNdim; k++) {
530 for (Int_t m=1; m<=fNdim; m++) {
531 if (fAdet->GetBinContent(k,m)) {
533 }
534 }
535 }
536
537 unfres = Unfold(GetKReg());
538
539 for (Int_t j=1; j<=fNdim; j++) {
540 toymean->SetBinContent(j, toymean->GetBinContent(j) + unfres->GetBinContent(j)/ntoys);
541 }
542 delete unfres;
543 unfres = nullptr;
544 }
545
546 // Reset the random seed
547 random.SetSeed(seed);
548
549 for (int i=1; i<=ntoys; i++) {
550 for (Int_t k=1; k<=fNdim; k++) {
551 for (Int_t m=1; m<=fNdim; m++) {
552 if (fAdet->GetBinContent(k,m))
554 }
555 }
556
557 unfres = Unfold(GetKReg());
558
559 for (Int_t j=1; j<=fNdim; j++) {
560 for (Int_t k=1; k<=fNdim; k++) {
561 unfcov->SetBinContent(j,k,unfcov->GetBinContent(j,k) + ( (unfres->GetBinContent(j) - toymean->GetBinContent(j))*(unfres->GetBinContent(k) - toymean->GetBinContent(k))/(ntoys-1)) );
562 }
563 }
564 delete unfres;
565 unfres = nullptr;
566 }
567 delete toymean;
569
570 return unfcov;
571}
572
573////////////////////////////////////////////////////////////////////////////////
574/// Returns d vector (for choosing appropriate regularisation)
575
577{
578 for (int i=1; i<=fDHist->GetNbinsX(); i++) {
580 }
581 return fDHist;
582}
583
584////////////////////////////////////////////////////////////////////////////////
585/// Returns singular values vector
586
588{
589 return fSVHist;
590}
591
592////////////////////////////////////////////////////////////////////////////////
593/// Returns the computed regularized covariance matrix corresponding to total uncertainties on measured spectrum as passed in the constructor.
594/// Note that this covariance matrix will not contain the effects of forced normalization if spectrum is normalized to unit area.
595
597{
598 return fXtau;
599}
600
601////////////////////////////////////////////////////////////////////////////////
602/// Returns the computed inverse of the covariance matrix
603
605{
606 return fXinv;
607}
608
609////////////////////////////////////////////////////////////////////////////////
610/// Returns the covariance matrix
611
613{
614 return fBcov;
615}
616
617////////////////////////////////////////////////////////////////////////////////
618/// Fill 1D histogram into vector
619
620void TSVDUnfold::H2V( const TH1D* histo, TVectorD& vec )
621{
622 for (Int_t i=0; i<histo->GetNbinsX(); i++) vec(i) = histo->GetBinContent(i+1);
623}
624
625////////////////////////////////////////////////////////////////////////////////
626/// Fill 1D histogram errors into vector
627
628void TSVDUnfold::H2Verr( const TH1D* histo, TVectorD& vec )
629{
630 for (Int_t i=0; i<histo->GetNbinsX(); i++) vec(i) = histo->GetBinError(i+1);
631}
632
633////////////////////////////////////////////////////////////////////////////////
634/// Fill vector into 1D histogram
635
636void TSVDUnfold::V2H( const TVectorD& vec, TH1D& histo )
637{
638 for(Int_t i=0; i<vec.GetNrows(); i++) histo.SetBinContent(i+1, vec(i));
639}
640
641////////////////////////////////////////////////////////////////////////////////
642/// Fill 2D histogram into matrix
643
644void TSVDUnfold::H2M( const TH2D* histo, TMatrixD& mat )
645{
646 for (Int_t j=0; j<histo->GetNbinsX(); j++) {
647 for (Int_t i=0; i<histo->GetNbinsY(); i++) {
648 mat(i,j) = histo->GetBinContent(i+1,j+1);
649 }
650 }
651}
652
653////////////////////////////////////////////////////////////////////////////////
654/// Fill 2D histogram into matrix
655
656void TSVDUnfold::M2H( const TMatrixD& mat, TH2D& histo )
657{
658 for (Int_t j=0; j<mat.GetNcols(); j++) {
659 for (Int_t i=0; i<mat.GetNrows(); i++) {
660 histo.SetBinContent(i+1,j+1, mat(i,j));
661 }
662 }
663}
664
665////////////////////////////////////////////////////////////////////////////////
666/// Divide entries of two vectors
667
668TVectorD TSVDUnfold::VecDiv( const TVectorD& vec1, const TVectorD& vec2, Int_t zero )
669{
670 TVectorD quot(vec1.GetNrows());
671 for (Int_t i=0; i<vec1.GetNrows(); i++) {
672 if (vec2(i) != 0) quot(i) = vec1(i) / vec2(i);
673 else {
674 if (zero) quot(i) = 0;
675 else quot(i) = vec1(i);
676 }
677 }
678 return quot;
679}
680
681////////////////////////////////////////////////////////////////////////////////
682/// Divide matrix entries by vector
683
685{
686 TMatrixD quotmat(mat.GetNrows(), mat.GetNcols());
687 for (Int_t i=0; i<mat.GetNrows(); i++) {
688 for (Int_t j=0; j<mat.GetNcols(); j++) {
689 if (vec(i) != 0) quotmat(i,j) = mat(i,j) / vec(i);
690 else {
691 if (zero) quotmat(i,j) = 0;
692 else quotmat(i,j) = mat(i,j);
693 }
694 }
695 }
696 return quotmat;
697}
698
699////////////////////////////////////////////////////////////////////////////////
700/// Multiply entries of two vectors
701
703{
704 TVectorD res(vec1.GetNrows());
705 for (Int_t i=0; i<vec1.GetNrows(); i++) res(i) = vec1(i) * vec2(i);
706 return res;
707}
708
709////////////////////////////////////////////////////////////////////////////////
710/// Compute curvature of vector
711
713{
714 return vec*(curv*vec);
715}
716
717////////////////////////////////////////////////////////////////////////////////
718
720{
721 Double_t eps = 0.00001;
722
723 Int_t ndim = tCurv.GetNrows();
724
725 // Init
726 tCurv *= 0;
727 tC *= 0;
728
729 if (fDdim == 0) for (Int_t i=0; i<ndim; i++) tC(i,i) = 1;
730 else if (fDdim == 1) {
731 for (Int_t i=0; i<ndim; i++) {
732 if (i < ndim-1) tC(i,i+1) = 1.0;
733 tC(i,i) = 1.0;
734 }
735 }
736 else if (fDdim == 2) {
737 for (Int_t i=0; i<ndim; i++) {
738 if (i > 0) tC(i,i-1) = 1.0;
739 if (i < ndim-1) tC(i,i+1) = 1.0;
740 tC(i,i) = -2.0;
741 }
742 tC(0,0) = -1.0;
743 tC(ndim-1,ndim-1) = -1.0;
744 }
745 else if (fDdim == 3) {
746 for (Int_t i=1; i<ndim-2; i++) {
747 tC(i,i-1) = 1.0;
748 tC(i,i) = -3.0;
749 tC(i,i+1) = 3.0;
750 tC(i,i+2) = -1.0;
751 }
752 }
753 else if (fDdim==4) {
754 for (Int_t i=0; i<ndim; i++) {
755 if (i > 0) tC(i,i-1) = -4.0;
756 if (i < ndim-1) tC(i,i+1) = -4.0;
757 if (i > 1) tC(i,i-2) = 1.0;
758 if (i < ndim-2) tC(i,i+2) = 1.0;
759 tC(i,i) = 6.0;
760 }
761 tC(0,0) = 2.0;
762 tC(ndim-1,ndim-1) = 2.0;
763 tC(0,1) = -3.0;
764 tC(ndim-2,ndim-1) = -3.0;
765 tC(1,0) = -3.0;
766 tC(ndim-1,ndim-2) = -3.0;
767 tC(1,1) = 6.0;
768 tC(ndim-2,ndim-2) = 6.0;
769 }
770 else if (fDdim == 5) {
771 for (Int_t i=2; i < ndim-3; i++) {
772 tC(i,i-2) = 1.0;
773 tC(i,i-1) = -5.0;
774 tC(i,i) = 10.0;
775 tC(i,i+1) = -10.0;
776 tC(i,i+2) = 5.0;
777 tC(i,i+3) = -1.0;
778 }
779 }
780 else if (fDdim == 6) {
781 for (Int_t i = 3; i < ndim - 3; i++) {
782 tC(i,i-3) = 1.0;
783 tC(i,i-2) = -6.0;
784 tC(i,i-1) = 15.0;
785 tC(i,i) = -20.0;
786 tC(i,i+1) = 15.0;
787 tC(i,i+2) = -6.0;
788 tC(i,i+3) = 1.0;
789 }
790 }
791
792 // Add epsilon to avoid singularities
793 for (Int_t i=0; i<ndim; i++) tC(i,i) = tC(i,i) + eps;
794
795 //Get curvature matrix
796 for (Int_t i=0; i<ndim; i++) {
797 for (Int_t j=0; j<ndim; j++) {
798 for (Int_t k=0; k<ndim; k++) {
799 tCurv(i,j) = tCurv(i,j) + tC(k,i)*tC(k,j);
800 }
801 }
802 }
803}
804
805////////////////////////////////////////////////////////////////////////////////
806
808{
809 fDHist = new TH1D( "dd", "d vector after orthogonal transformation", fNdim, 0, fNdim );
810 fDHist->Sumw2();
811
812 fSVHist = new TH1D( "sv", "Singular values of AC^-1", fNdim, 0, fNdim );
813 fSVHist->Sumw2();
814
815 fXtau = (TH2D*)fAdet->Clone("Xtau");
816 fXtau->SetTitle("Regularized covariance matrix");
817 fXtau->Sumw2();
818
819 fXinv = (TH2D*)fAdet->Clone("Xinv");
820 fXinv->SetTitle("Inverse covariance matrix");
821 fXinv->Sumw2();
822}
823
824////////////////////////////////////////////////////////////////////////////////
825/// naive regularised inversion cuts off small elements
826
828{
829 // init reduced matrix
830 const UInt_t n = mat.GetNrows();
831 UInt_t nn = 0;
832
833 UInt_t *ipos = new UInt_t[n];
834 // UInt_t ipos[n];
835
836 // find max diagonal entries
837 Double_t ymax = 0;
838 for (UInt_t i=0; i<n; i++) if (TMath::Abs(mat[i][i]) > ymax) ymax = TMath::Abs(mat[i][i]);
839
840 for (UInt_t i=0; i<n; i++) {
841
842 // save position of accepted entries
843 if (TMath::Abs(mat[i][i])/ymax > eps) ipos[nn++] = i;
844 }
845
846 // effective matrix
847 TMatrixDSym matwork( nn );
848 for (UInt_t in=0; in<nn; in++) for (UInt_t jn=0; jn<nn; jn++) matwork(in,jn) = 0;
849
850 // fill non-zero effective working matrix
851 for (UInt_t in=0; in<nn; in++) {
852
853 matwork[in][in] = mat[ipos[in]][ipos[in]];
854 for (UInt_t jn=in+1; jn<nn; jn++) {
855 matwork[in][jn] = mat[ipos[in]][ipos[jn]];
856 matwork[jn][in] = matwork[in][jn];
857 }
858 }
859
860 // invert
861 matwork.Invert();
862
863 // reinitialise old matrix
864 for (UInt_t i=0; i<n; i++) for (UInt_t j=0; j<n; j++) mat[i][j] = 0;
865
866 // refill inverted matrix in old one
867 for (UInt_t in=0; in<nn; in++) {
868 mat[ipos[in]][ipos[in]] = matwork[in][in];
869 for (UInt_t jn=in+1; jn<nn; jn++) {
870 mat[ipos[in]][ipos[jn]] = matwork[in][jn];
871 mat[ipos[jn]][ipos[in]] = mat[ipos[in]][ipos[jn]];
872 }
873 }
874 delete [] ipos;
875}
876
877////////////////////////////////////////////////////////////////////////////////
878/// Helper routine to compute chi-squared between distributions using the computed inverse of the covariance matrix for the unfolded spectrum as given in paper.
879
880Double_t TSVDUnfold::ComputeChiSquared( const TH1D& truspec, const TH1D& unfspec)
881{
882 UInt_t n = truspec.GetNbinsX();
883
884 // compute chi2
885 Double_t chi2 = 0;
886 for (UInt_t i=0; i<n; i++) {
887 for (UInt_t j=0; j<n; j++) {
888 chi2 += ( (truspec.GetBinContent( i+1 )-unfspec.GetBinContent( i+1 )) *
889 (truspec.GetBinContent( j+1 )-unfspec.GetBinContent( j+1 )) * fXinv->GetBinContent(i+1,j+1) );
890 }
891 }
892
893 return chi2;
894}
895
#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:94
#define ClassImp(name)
Definition Rtypes.h:382
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:670
void SetTitle(const char *title) override
Change/set the title.
Definition TH1.cxx:6718
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:9063
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:9206
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:9222
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition TH1.cxx:5061
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:9020
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:1021
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:967
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