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