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