Logo ROOT  
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  \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 using 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 
77 TSVDUnfold::TSVDUnfold( const TH1D *bdat, const TH1D *bini, const TH1D *xini, const TH2D *Adet )
78 : TObject (),
79 fNdim (0),
80 fDdim (2),
81 fNormalize (kFALSE),
82 fKReg (-1),
83 fDHist (NULL),
84 fSVHist (NULL),
85 fXtau (NULL),
86 fXinv (NULL),
87 fBdat (bdat),
88 fBini (bini),
89 fXini (xini),
90 fAdet (Adet),
91 fToyhisto (NULL),
92 fToymat (NULL),
93 fToyMode (kFALSE),
94 fMatToyMode (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 
130 TSVDUnfold::TSVDUnfold( const TH1D *bdat, TH2D* Bcov, const TH1D *bini, const TH1D *xini, const TH2D *Adet )
131 : TObject (),
132 fNdim (0),
133 fDdim (2),
134 fNormalize (kFALSE),
135 fKReg (-1),
136 fDHist (NULL),
137 fSVHist (NULL),
138 fXtau (NULL),
139 fXinv (NULL),
140 fBdat (bdat),
141 fBcov (Bcov),
142 fBini (bini),
143 fXini (xini),
144 fAdet (Adet),
145 fToyhisto (NULL),
146 fToymat (NULL),
147 fToyMode (kFALSE),
148 fMatToyMode (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 ),
177 fNdim (other.fNdim),
178 fDdim (other.fDdim),
179 fNormalize (other.fNormalize),
180 fKReg (other.fKReg),
181 fDHist (other.fDHist),
182 fSVHist (other.fSVHist),
183 fXtau (other.fXtau),
184 fXinv (other.fXinv),
185 fBdat (other.fBdat),
186 fBcov (other.fBcov),
187 fBini (other.fBini),
188 fXini (other.fXini),
189 fAdet (other.fAdet),
190 fToyhisto (other.fToyhisto),
191 fToymat (other.fToymat),
192 fToyMode (other.fToyMode),
193 fMatToyMode (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 histogams 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
453  TVectorD g(fNdim);
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
481  TVectorD g(fNdim);
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++) {
489  fToyhisto->SetBinContent( j, fBdat->GetBinContent(j)+g(j-1) );
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;
504  fToyMode = kFALSE;
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)) {
536  fToymat->SetBinContent(k, m, random.Poisson(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))
557  fToymat->SetBinContent(k, m, random.Poisson(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 
624 void 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 
632 void 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 
640 void 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 
648 void 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 
660 void 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 
672 TVectorD 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 
688 TMatrixD TSVDUnfold::MatDivVec( const TMatrixD& mat, const TVectorD& vec, Int_t zero )
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 
706 TVectorD TSVDUnfold::CompProd( const TVectorD& vec1, const TVectorD& vec2 )
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 
884 Double_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 
TVectorT::GetNrows
Int_t GetNrows() const
Definition: TVectorT.h:81
TDecompSVD::GetU
const TMatrixD & GetU()
Definition: TDecompSVD.h:61
TSVDUnfold::fDHist
TH1D * fDHist
Regularisation parameter.
Definition: TSVDUnfold.h:131
m
auto * m
Definition: textangle.C:8
n
const Int_t n
Definition: legend1.C:16
TSVDUnfold::InitHistos
void InitHistos()
Definition: TSVDUnfold.cxx:811
TSVDUnfold::ComputeChiSquared
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:884
TRandom::Gaus
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:263
ymax
float ymax
Definition: THbookFile.cxx:95
e
#define e(i)
Definition: RSha256.hxx:121
TSVDUnfold::GetD
TH1D * GetD() const
Returns d vector (for choosing appropriate regularisation)
Definition: TSVDUnfold.cxx:580
TSVDUnfold::GetSV
TH1D * GetSV() const
Returns singular values vector.
Definition: TSVDUnfold.cxx:591
TSVDUnfold::M2H
static void M2H(const TMatrixD &mat, TH2D &histo)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:660
TDecompSVD
Definition: TDecompSVD.h:23
TH2::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:88
TSVDUnfold::fToymat
TH2D * fToymat
Toy MC histogram.
Definition: TSVDUnfold.h:145
ClassImp
#define ClassImp(name)
Definition: Rtypes.h:364
Form
char * Form(const char *fmt,...)
TSVDUnfold::RegularisedSymMatInvert
static void RegularisedSymMatInvert(TMatrixDSym &mat, Double_t eps=1e-3)
naive regularised inversion cuts off small elements
Definition: TSVDUnfold.cxx:831
TSVDUnfold::V2H
static void V2H(const TVectorD &vec, TH1D &histo)
Fill vector into 1D histogram.
Definition: TSVDUnfold.cxx:640
TObject::Info
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:864
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:691
TH1D
1-D histogram with a double per channel (see TH1 documentation)}
Definition: TH1.h:615
TSVDUnfold.h
TSVDUnfold::fBdat
const TH1D * fBdat
Computed inverse of covariance matrix.
Definition: TSVDUnfold.h:137
TObject::Fatal
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:918
TH1::GetBinError
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8518
TH1::SetBinContent
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:8677
TSVDUnfold::fMatToyMode
Bool_t fMatToyMode
Internal switch for covariance matrix propagation.
Definition: TSVDUnfold.h:147
TSVDUnfold::H2M
static void H2M(const TH2D *histo, TMatrixD &mat)
Fill 2D histogram into matrix.
Definition: TSVDUnfold.cxx:648
TMatrixTSym
Definition: TMatrixDSymfwd.h:22
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120
TSVDUnfold::GetAdetCovMatrix
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:515
TString
Definition: TString.h:136
TDecompSVD.h
TMatrixT
Definition: TMatrixDfwd.h:22
TDecompSVD::GetSig
const TVectorD & GetSig()
Definition: TDecompSVD.h:65
TSVDUnfold::fXtau
TH2D * fXtau
Distribution of singular values.
Definition: TSVDUnfold.h:133
TSVDUnfold::fXini
const TH1D * fXini
Definition: TSVDUnfold.h:140
TH1::SetTitle
virtual void SetTitle(const char *title)
See GetStatOverflows for more information.
Definition: TH1.cxx:6344
TH1::Clone
TObject * Clone(const char *newname=0) const
Make a complete copy of the underlying object.
Definition: TH1.cxx:2664
TGeant4Unit::L
static constexpr double L
Definition: TGeant4SystemOfUnits.h:123
TSVDUnfold::GetCurvature
static Double_t GetCurvature(const TVectorD &vec, const TMatrixD &curv)
Compute curvature of vector.
Definition: TSVDUnfold.cxx:716
TSVDUnfold::CompProd
static TVectorD CompProd(const TVectorD &vec1, const TVectorD &vec2)
Multiply entries of two vectors.
Definition: TSVDUnfold.cxx:706
TDecompSVD::GetV
const TMatrixD & GetV()
Definition: TDecompSVD.h:63
TSVDUnfold::fKReg
Int_t fKReg
Normalize unfolded spectrum to 1.
Definition: TSVDUnfold.h:130
TSVDUnfold::fToyhisto
TH1D * fToyhisto
Definition: TSVDUnfold.h:144
TH1::GetBinContent
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4906
TMatrixT::kTransposed
@ kTransposed
Definition: TMatrixT.h:58
TSVDUnfold::GetXinv
TH2D * GetXinv() const
Returns the computed inverse of the covariance matrix.
Definition: TSVDUnfold.cxx:608
TH2D.h
TRandom3
Definition: TRandom3.h:27
h
#define h(i)
Definition: RSha256.hxx:124
TRandom3::SetSeed
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:206
a
auto * a
Definition: textangle.C:12
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:92
TH2D
2-D histogram with a double per channel (see TH1 documentation)}
Definition: TH2.h:292
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:735
TRandom3.h
TSVDUnfold::fAdet
const TH2D * fAdet
Definition: TSVDUnfold.h:141
TSVDUnfold::GetBCov
TH2D * GetBCov() const
Returns the covariance matrix.
Definition: TSVDUnfold.cxx:616
TSVDUnfold::fXinv
TH2D * fXinv
Computed regularized covariance matrix.
Definition: TSVDUnfold.h:134
TSVDUnfold::fDdim
Int_t fDdim
Truth and reconstructed dimensions.
Definition: TSVDUnfold.h:128
unsigned int
TSVDUnfold::GetXtau
TH2D * GetXtau() const
Returns the computed regularized covariance matrix corresponding to total uncertainties on measured s...
Definition: TSVDUnfold.cxx:600
TSVDUnfold
Definition: TSVDUnfold.h:46
TSVDUnfold::~TSVDUnfold
virtual ~TSVDUnfold()
Destructor.
Definition: TSVDUnfold.cxx:200
TSVDUnfold::fToyMode
Bool_t fToyMode
Toy MC detector response matrix.
Definition: TSVDUnfold.h:146
TSVDUnfold::H2V
static void H2V(const TH1D *histo, TVectorD &vec)
Fill 1D histogram into vector.
Definition: TSVDUnfold.cxx:624
TVectorT
Definition: TMatrixTBase.h:78
Double_t
double Double_t
Definition: RtypesCore.h:59
TVectorT::Sum
Element Sum() const
Compute sum of elements.
Definition: TVectorT.cxx:637
TMatrixTSym::Invert
TMatrixTSym< Element > & Invert(Double_t *det=0)
Invert the matrix and calculate its determinant Notice that the LU decomposition is used instead of B...
Definition: TMatrixTSym.cxx:961
TH1::SetBinError
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:8661
TMatrixD
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
TRandom::Poisson
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
Definition: TRandom.cxx:391
TH1::Sumw2
virtual void Sumw2(Bool_t flag=kTRUE)
Create structure to store sum of squares of weights.
Definition: TH1.cxx:8475
TMatrixT::Transpose
TMatrixT< Element > & Transpose(const TMatrixT< Element > &source)
Transpose matrix source.
Definition: TMatrixT.cxx:1470
TObject
Definition: TObject.h:37
TSVDUnfold::H2Verr
static void H2Verr(const TH1D *histo, TVectorD &vec)
Fill 1D histogram errors into vector.
Definition: TSVDUnfold.cxx:632
TSVDUnfold::fSVHist
TH1D * fSVHist
Distribution of d (for checking regularization)
Definition: TSVDUnfold.h:132
TSVDUnfold::Unfold
TH1D * Unfold(Int_t kreg)
Perform the unfolding with regularisation parameter kreg.
Definition: TSVDUnfold.cxx:241
TSVDUnfold::GetUnfoldCovMatrix
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:409
TSVDUnfold::fNdim
Int_t fNdim
Definition: TSVDUnfold.h:127
TH1::GetNbinsY
virtual Int_t GetNbinsY() const
Definition: TH1.h:294
TSVDUnfold::MatDivVec
static TMatrixD MatDivVec(const TMatrixD &mat, const TVectorD &vec, Int_t zero=0)
Divide matrix entries by vector.
Definition: TSVDUnfold.cxx:688
TSVDUnfold::fBcov
TH2D * fBcov
Definition: TSVDUnfold.h:138
TSVDUnfold::TSVDUnfold
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:77
TH1D.h
TSVDUnfold::FillCurvatureMatrix
void FillCurvatureMatrix(TMatrixD &tCurv, TMatrixD &tC) const
Definition: TSVDUnfold.cxx:723
TSVDUnfold::VecDiv
static TVectorD VecDiv(const TVectorD &vec1, const TVectorD &vec2, Int_t zero=0)
Divide entries of two vectors.
Definition: TSVDUnfold.cxx:672
TSVDUnfold::fBini
const TH1D * fBini
Definition: TSVDUnfold.h:139
TSVDUnfold::GetKReg
Int_t GetKReg() const
Definition: TSVDUnfold.h:86
TH2::SetBinContent
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2480
TMath.h
TH1::GetNbinsX
virtual Int_t GetNbinsX() const
Definition: TH1.h:293
int
TSVDUnfold::fNormalize
Bool_t fNormalize
Derivative for curvature matrix.
Definition: TSVDUnfold.h:129
g
#define g(i)
Definition: RSha256.hxx:123