Logo ROOT  
Reference Guide
 
Loading...
Searching...
No Matches
TUnfoldSys.cxx
Go to the documentation of this file.
1// Author: Stefan Schmitt
2// DESY, 23/01/09
3
4// Version 17.9, parallel to changes in TUnfold
5//
6// History:
7// Version 17.8, parallel to changes in TUnfold
8// Version 17.7, bug fix in GetBackground()
9// Version 17.6, with updated doxygen comments
10// Version 17.5, in parallel to changes in TUnfold
11// Version 17.4, in parallel to changes in TUnfoldBinning
12// Version 17.3, in parallel to changes in TUnfoldBinning
13// Version 17.2, add methods to find back systematic and background sources
14// Version 17.1, bug fix with background uncertainty
15// Version 17.0, possibility to specify an error matrix with SetInput
16// Version 16.1, parallel to changes in TUnfold
17// Version 16.0, parallel to changes in TUnfold
18// Version 15, fix bugs with uncorr. uncertainties, add backgnd subtraction
19// Version 14, remove some print-out, do not add unused sys.errors
20// Version 13, support for systematic errors
21
22/** \class TUnfoldSys
23\ingroup Unfold
24An algorithm to unfold distributions from detector to truth level,
25with background subtraction and propagation of systematic uncertainties
26
27TUnfoldSys is used to decompose a measurement y into several sources x,
28given the measurement uncertainties, background b and a matrix of migrations A.
29The method can be applied to a large number of problems,
30where the measured distribution y is a linear superposition
31of several Monte Carlo shapes. Beyond such a simple template fit,
32TUnfoldSys has an adjustable regularisation term and also supports an
33optional constraint on the total number of events.
34Background sources can be specified, with a normalisation constant and
35normalisation uncertainty. In addition, variants of the response
36matrix may be specified, these are taken to determine systematic
37uncertainties.
38
39<b>For most applications, it is better to use the derived class
40TUnfoldDensity instead of TUnfoldSys. TUnfoldDensity adds
41features to TUnfoldSys, related to possible complex multidimensional
42arrangements of bins. For innocent
43users, the most notable improvement of TUnfoldDensity over TUnfoldSys are
44the getter functions. For TUnfoldSys, histograms have to be booked by the
45user and the getter functions fill the histogram bins. TUnfoldDensity
46simply returns a new, already filled histogram.</b>
47
48If you use this software, please consider the following citation
49
50<b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
51
52Detailed documentation and updates are available on
53http://www.desy.de/~sschmitt
54
55Brief recipy to use TUnfoldSys:
56<ul>
57<li>a matrix (truth,reconstructed) is given as a two-dimensional histogram
58 as argument to the constructor of TUnfold</li>
59<li>a vector of measurements is given as one-dimensional histogram using
60 the SetInput() method</li>
61<li>repeated calls to SubtractBackground() to specify background
62sources</li>
63<li>repeated calls to AddSysError() to specify systematic uncertainties
64<li>The unfolding is performed
65<ul>
66<li>either once with a fixed parameter tau, method DoUnfold(tau)</li>
67<li>or multiple times in a scan to determine the best chouce of tau,
68 method ScanLCurve()</li>
69</ul>
70<li>Unfolding results are retrieved using various GetXXX() methods
71</ul>
72
73Description of (systematic) uncertainties available in
74TUnfoldSys. There are covariance matrix contributions and there are
75systematic shifts. Systematic shifts correspond to the variation of a
76(buicance) parameter, for example a background normalisation or a
77one-sigma variation of a correlated systematic error.
78<table>
79<tr><th> </th><th>Set by</th>
80<th>Access covariance matrix</th>
81<th>Access vector of shifts</th>
82<th>Description</th>
83</tr>
84<tr>
85<td>(a)</td><td>TUnfoldSys constructor</td>
86<td>GetEmatrixSysUncorr()</td><td> n.a. </td>
87<td>uncorrelated errors on the input matrix histA, taken as the errors
88provided with the histogram. These are typically statistical errors
89from finite Monte Carlo samples.</td>
90</tr>
91<tr>
92<td>(b)</td><td>AddSysError()</td><td>GetEmatrixSysSource()</td>
93<td>GetDeltaSysSource()</td>
94<td>correlated shifts of the input matrix histA. These shifts are taken
95as one-sigma effects when switchig on a given error soure. Several
96such error sources may be defined</td>
97</tr>
98<tr>
99<td>(c)</td><td>SetTauError()</td><td>GetEmatrixSysTau()</td>
100<td>GetDeltaSysTau()</td>
101<td>A systematic error on the regularisation parameter tau</td>
102</tr>
103<tr><td>(d)</td><td>SubtractBackground()</td>
104<td>GetEmatrixSysBackgroundUncorr()</td><td>n.a.</td>
105<td>uncorrelated errors on background sources, originating from the errors
106 provided with the background histograms</td>
107</tr>
108<tr>
109<td>(e)</td><td>SubtractBackground</td>
110<td>GetEmatrixSysBackgroundScale()</td><td>GetDeltaSysBackgroundScale()</td>
111<td>scale errors on background sources</td>
112</tr>
113<tr>
114<td>(i)</td><td>SetInput()</td>
115<td>GetEmatrixInput()</td><td>n.a.</td><td>statistical uncertainty of
116the input (the measurement)</td>
117</tr>
118<tr>
119<td>(i)+(d)+(e)</td><td>see above</td><td>GetEmatrix()</td><td>n.a.</td>
120<td>Partial sun of uncertainties: all sources which are propagated to the covariance before unfolding</td>
121</tr>
122<tr>
123<td>
124(i)+(a)+(b)+(c)+(d)+(e)</td><td>see above</td><td>GetEmatrixTotal()</td>
125<td>n.a.</td><td>All known error sources summed up</td>
126</tr>
127</table>
128
129Note: (a), (b), (c) are propagated to the result AFTER unfolding,
130whereas the background errors (d) and (e) are added to the data errors
131BEFORE unfolding. For this reason the errors of type (d) and (e) are
132INCLUDED in the standard error matrix and other methods provided by
133the base class TUnfold, whereas errors of type (a), (b), (c) are NOT
134INCLUDED in the methods provided by the base class TUnfold.
135*/
136
137/*
138 This file is part of TUnfold.
139
140 TUnfold is free software: you can redistribute it and/or modify
141 it under the terms of the GNU General Public License as published by
142 the Free Software Foundation, either version 3 of the License, or
143 (at your option) any later version.
144
145 TUnfold is distributed in the hope that it will be useful,
146 but WITHOUT ANY WARRANTY; without even the implied warranty of
147 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
148 GNU General Public License for more details.
149
150 You should have received a copy of the GNU General Public License
151 along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
152*/
153
154#include <iostream>
155#include <TMap.h>
156#include <TMath.h>
157#include <TObjString.h>
158#include <TSortedList.h>
159#include <RVersion.h>
160#include <cmath>
161
162#include "TUnfoldSys.h"
163
164
166{
167 // delete all data members
170 delete fBgrIn;
171 delete fBgrErrUncorrInSq;
172 delete fBgrErrScaleIn;
173 delete fSysIn;
174 ClearResults();
175 delete fDeltaCorrX;
176 delete fDeltaCorrAx;
179}
180
181////////////////////////////////////////////////////////////////////////
182/// only for use by root streamer or derived classes
183///
185{
186 // set all pointers to zero
188}
189
190////////////////////////////////////////////////////////////////////////
191/// set up response matrix A, uncorrelated uncertainties of A and
192/// regularisation scheme
193///
194/// \param[in] hist_A matrix that describes the migrations
195/// \param[in] histmap mapping of the histogram axes to the unfolding output
196/// \param[in] regmode (default=kRegModeSize) global regularisation mode
197/// \param[in] constraint (default=kEConstraintArea) type of constraint
198///
199/// For further details, consult the constructir of the class TUnfold.
200/// The uncertainties of hist_A are taken to be uncorrelated and aper
201/// propagated to the unfolding result, method GetEmatrixSysUncorr().
204 : TUnfold(hist_A,histmap,regmode,constraint)
205{
206 // data members initialized to something different from zero:
207 // fDA2, fDAcol
208
209 // initialize TUnfoldSys
211
212 // save underflow and overflow bins
213 fAoutside = new TMatrixD(GetNx(),2);
214 // save the normalized errors on hist_A
215 // to the matrices fDAinRelSq and fDAinColRelSq
216 fDAinColRelSq = new TMatrixD(GetNx(),1);
217
218 Int_t nmax=GetNx()*GetNy();
222
224 for (Int_t ix = 0; ix < GetNx(); ix++) {
225 Int_t ibinx = fXToHist[ix];
227 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
228 Double_t dz;
230 dz = hist_A->GetBinError(ibinx, ibiny);
231 } else {
232 dz = hist_A->GetBinError(ibiny, ibinx);
233 }
235 // quadratic sum of all errors from all bins,
236 // including under/overflow bins
237 (*fDAinColRelSq)(ix,0) += normerr_sq;
238
239 if(ibiny==0) {
240 // underflow bin
242 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibinx, ibiny);
243 } else {
244 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibiny, ibinx);
245 }
246 } else if(ibiny==GetNy()+1) {
247 // overflow bins
249 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibinx, ibiny);
250 } else {
251 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibiny, ibinx);
252 }
253 } else {
254 // error on this bin
259 }
260 }
261 }
262 if(da_nonzero) {
265 } else {
267 }
268 delete[] rowDAinRelSq;
269 delete[] colDAinRelSq;
270 delete[] dataDAinRelSq;
271}
272
273////////////////////////////////////////////////////////////////////////
274/// Specify a correlated systematic uncertainty
275///
276/// \param[in] sysError alternative matrix or matrix of absolute/relative shifts
277/// \param[in] name identifier of the error source
278/// \param[in] histmap mapping of the histogram axes
279/// \param[in] mode format of the error source
280///
281/// <b>sysError</b> corresponds to a one-sigma variation. If
282/// may be given in various forms, specified by <b>mode</b>
283/// <ul>
284/// <li><b>mode=kSysErrModeMatrix</b> the histogram <b>sysError</b>
285/// corresponds to an alternative response matrix.</li>
286/// <li><b>mode=kSysErrModeShift</b> the content of the histogram <b>sysError</b> are the absolute shifts of the response matrix
287/// <li><b>mode=kSysErrModeRelative</b> the content of the histogram <b>sysError</b>
288/// specifies the relative uncertainties
289/// </ul>
290/// Internally, all three cases are transformed to the case <b>mode=kSysErrModeMatrix</b>.
292(const TH2 *sysError,const char *name,EHistMap histmap,ESysErrMode mode)
293{
294
295 if(fSysIn->FindObject(name)) {
296 Error("AddSysError","Source %s given twice, ignoring 2nd call.\n",name);
297 } else {
298 // a copy of fA is made. It can be accessed inside the loop
299 // without having to take care that the sparse structure of *fA
300 // otherwise, *fA may be accidentally destroyed by asking
301 // for an element which is zero.
302 TMatrixD aCopy(*fA);
303
304 Int_t nmax= GetNx()*GetNy();
306 Int_t *cols=new Int_t[nmax];
307 Int_t *rows=new Int_t[nmax];
308 nmax=0;
309 for (Int_t ix = 0; ix < GetNx(); ix++) {
310 Int_t ibinx = fXToHist[ix];
311 Double_t sum=0.0;
312 for(Int_t loop=0;loop<2;loop++) {
313 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
314 Double_t z;
315 // get bin content, depending on histmap
317 z = sysError->GetBinContent(ibinx, ibiny);
318 } else {
319 z = sysError->GetBinContent(ibiny, ibinx);
320 }
321 // correct to absolute numbers
323 Double_t z0;
324 if((ibiny>0)&&(ibiny<=GetNy())) {
325 z0=aCopy(ibiny-1,ix)*fSumOverY[ix];
326 } else if(ibiny==0) {
327 z0=(*fAoutside)(ix,0);
328 } else {
329 z0=(*fAoutside)(ix,1);
330 }
332 z += z0;
333 } else if(mode==kSysErrModeRelative) {
334 z = z0*(1.+z);
335 }
336 }
337 if(loop==0) {
338 // sum up all entries, including overflow bins
339 sum += z;
340 } else {
341 if((ibiny>0)&&(ibiny<=GetNy())) {
342 // save normalized matrix of shifts,
343 // excluding overflow bins
344 rows[nmax]=ibiny-1;
345 cols[nmax]=ix;
346 if(sum>0.0) {
347 data[nmax]=z/sum-aCopy(ibiny-1,ix);
348 } else {
349 data[nmax]=0.0;
350 }
351 if(data[nmax] != 0.0) nmax++;
352 }
353 }
354 }
355 }
356 }
357 if(nmax==0) {
358 Error("AddSysError",
359 "source %s has no influence and has not been added.\n",name);
360 } else {
364 }
365 delete[] data;
366 delete[] rows;
367 delete[] cols;
368 }
369}
370
371////////////////////////////////////////////////////////////////////////
372/// perform background subtraction
373///
374/// This prepares the data members for the base class TUnfold, such
375/// that the background is properly taken into account.
377{
378 // fY = fYData - fBgrIn
379 // fVyy = fVyyData + fBgrErrUncorr^2 + fBgrErrCorr * fBgrErrCorr#
380 // fVyyinv = fVyy^(-1)
381
382 if(fYData) {
385 if(fBgrIn->GetEntries()<=0) {
386 // simple copy
387 fY=new TMatrixD(*fYData);
389 } else {
390 if(GetVyyInv()) {
391 Warning("DoBackgroundSubtraction",
392 "inverse error matrix from user input,"
393 " not corrected for background");
394 }
395 // copy of the data
396 fY=new TMatrixD(*fYData);
397 // subtract background from fY
398 const TObject *key;
399 {
401 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
402 const TMatrixD *bgr=(const TMatrixD *)((const TPair *)*bgrPtr)->Value();
403 for(Int_t i=0;i<GetNy();i++) {
404 (*fY)(i,0) -= (*bgr)(i,0);
405 }
406 }
407 }
408 // copy original error matrix
410 // determine used bins
415 Int_t *usedBin=new Int_t[ny];
416 for(Int_t i=0;i<ny;i++) {
417 usedBin[i]=0;
418 }
419 for(Int_t i=0;i<ny;i++) {
420 for(Int_t k=vyydata_rows[i];k<vyydata_rows[i+1];k++) {
421 if(vyydata_data[k]>0.0) {
422 usedBin[i]++;
423 usedBin[vyydata_cols[k]]++;
424 }
425 }
426 }
427 // add uncorrelated background errors
428 {
430 for(key=bgrErrUncorrSqPtr.Next();key;
431 key=bgrErrUncorrSqPtr.Next()) {
432 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)((const TPair *)*bgrErrUncorrSqPtr)->Value();
433 for(Int_t yi=0;yi<ny;yi++) {
434 if(!usedBin[yi]) continue;
435 vyy(yi,yi) +=(*bgrerruncorrSquared)(yi,0);
436 }
437 }
438 }
439 // add correlated background errors
440 {
442 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
443 const TMatrixD *bgrerrscale=(const TMatrixD *)((const TPair *)*bgrErrScalePtr)->Value();
444 for(Int_t yi=0;yi<ny;yi++) {
445 if(!usedBin[yi]) continue;
446 for(Int_t yj=0;yj<ny;yj++) {
447 if(!usedBin[yj]) continue;
448 vyy(yi,yj) +=(*bgrerrscale)(yi,0)* (*bgrerrscale)(yj,0);
449 }
450 }
451 }
452 }
453 delete[] usedBin;
454 usedBin=nullptr;
455
456 // convert to sparse matrix
458
459 }
460 } else {
461 Fatal("DoBackgroundSubtraction","No input vector defined");
462 }
463}
464
467 const TH2 *hist_vyy_inv)
468{
469 // Define the input data for subsequent calls to DoUnfold(Double_t)
470 // input: input distribution with errors
471 // scaleBias: scale factor applied to the bias
472 // oneOverZeroError: for bins with zero error, this number defines 1/error.
473 // Return value: number of bins with bad error
474 // +10000*number of unconstrained output bins
475 // Note: return values>=10000 are fatal errors,
476 // for the given input, the unfolding can not be done!
477 // Calls the SetInput method of the base class, then renames the input
478 // vectors fY and fVyy, then performs the background subtraction
479 // Data members modified:
480 // fYData,fY,fVyyData,fVyy,fVyyinvData,fVyyinv
481 // and those modified by TUnfold::SetInput()
482 // and those modified by DoBackgroundSubtraction()
483
486 fYData=fY;
487 fY=nullptr;
489 fVyy=nullptr;
491
492 return r;
493}
494
495////////////////////////////////////////////////////////////////////////
496/// Specify a source of background
497///
498/// \param[in] bgr background distribution with uncorrelated errors
499/// \param[in] name identifier for this background source
500/// \param[in] scale normalisation factor applied to the background
501/// \param[in] scaleError normalisation uncertainty
502///
503/// The contribution <b>scale</b>*<b>bgr</b> is subtracted from the
504/// measurement prior to unfolding. The following contributions are
505/// added to the input covarianc ematrix
506/// <ul>
507/// <li>using the uncorrelated histogram errors <b>dbgr</b>, the contribution
508/// (<b>scale</b>*<b>dbgr<sub>i</sub></b>)<sup>2</sup> is added to the
509/// diagonals of the covariance</li>
510/// <li>using the histogram contents, the background normalisation uncertainty contribution
511/// <b>dscale</b>*<b>bgr<sub>i</sub></b> <b>dscale</b>*<b>bgr<sub>j</sub></b>
512/// is added to the covariance matrix</li>
513
515(const TH1 *bgr,const char *name,Double_t scale,Double_t scale_error)
516{
517 // Data members modified:
518 // fBgrIn,fBgrErrUncorrInSq,fBgrErrScaleIn
519 // and those modified by DoBackgroundSubtraction()
520
521 // save background source
522 if(fBgrIn->FindObject(name)) {
523 Error("SubtractBackground","Source %s given twice, ignoring 2nd call.\n",
524 name);
525 } else {
529 for(Int_t row=0;row<GetNy();row++) {
530 (*bgrScaled)(row,0) = scale*bgr->GetBinContent(row+1);
531 (*bgrErrUncSq)(row,0) =
532 TMath::Power(scale*bgr->GetBinError(row+1),2.);
533 (*bgrErrCorr)(row,0) = scale_error*bgr->GetBinContent(row+1);
534 }
538 if(fYData) {
540 } else {
541 Info("SubtractBackground",
542 "Background subtraction prior to setting input data");
543 }
544 }
545}
546
547////////////////////////////////////////////////////////////////////////
548/// get background into a histogram
549///
550/// \param[inout] bgrHist target histogram, content and errors will be altered
551/// \param[in] bgrSource (default=nullptr) name of backgrond source or zero
552/// to add all sources of background
553/// \param[in] binMap (default=nullptr) remap histogram bins
554/// \param[in] includeError (default=3) include uncorrelated(1),
555/// correlated (2) or both (3) sources of uncertainty in the
556/// histogram errors
557/// \param[in] clearHist (default=true) reset histogram before adding
558/// up the specified background sources
559///
560/// the array <b>binMap</b> is explained with the method GetOutput().
561/// The flag <b>clearHist</b> may be used to add background from
562/// several sources in successive calls to GetBackground().
563
565(TH1 *bgrHist,const char *bgrSource,const Int_t *binMap,
567{
568 if(clearHist) {
570 }
571 // add all background sources
572 const TObject *key;
573 {
575 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
576 TString bgrName=((const TObjString *)key)->GetString();
577 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
578 const TMatrixD *bgr=(const TMatrixD *)((const TPair *)*bgrPtr)->Value();
579 for(Int_t i=0;i<GetNy();i++) {
580 Int_t destBin=binMap[i+1];
581 bgrHist->SetBinContent(destBin,bgrHist->GetBinContent(destBin)+
582 (*bgr)(i,0));
583 }
584 }
585 }
586 // add uncorrelated background errors
587 if(includeError &1) {
589 for(key=bgrErrUncorrSqPtr.Next();key;key=bgrErrUncorrSqPtr.Next()) {
590 TString bgrName=((const TObjString *)key)->GetString();
591 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
592 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)((const TPair *)*bgrErrUncorrSqPtr)->Value();
593 for(Int_t i=0;i<GetNy();i++) {
594 Int_t destBin=binMap[i+1];
595 bgrHist->SetBinError
597 ((*bgrerruncorrSquared)(i,0)+
598 TMath::Power(bgrHist->GetBinError(destBin),2.)));
599 }
600 }
601 }
602 if(includeError & 2) {
604 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
605 TString bgrName=((const TObjString *)key)->GetString();
606 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
607 const TMatrixD *bgrerrscale=(TMatrixD const *)((const TPair *)*bgrErrScalePtr)->Value();
608 for(Int_t i=0;i<GetNy();i++) {
609 Int_t destBin=binMap[i+1];
610 bgrHist->SetBinError(destBin,hypot((*bgrerrscale)(i,0),
611 bgrHist->GetBinError(destBin)));
612 }
613 }
614 }
615}
616
618{
619 // initialize pointers and TMaps
620
621 // input
622 fDAinRelSq = nullptr;
623 fDAinColRelSq = nullptr;
624 fAoutside = nullptr;
625 fBgrIn = new TMap();
626 fBgrErrUncorrInSq = new TMap();
627 fBgrErrScaleIn = new TMap();
628 fSysIn = new TMap();
633 // results
634 fEmatUncorrX = nullptr;
635 fEmatUncorrAx = nullptr;
636 fDeltaCorrX = new TMap();
637 fDeltaCorrAx = new TMap();
640 fDeltaSysTau = nullptr;
641 fDtau=0.0;
642 fYData=nullptr;
643 fVyyData=nullptr;
644}
645
646////////////////////////////////////////////////////////////////////////////////
647/// Clear all data members which depend on the unfolding results.
648
650{
651 // clear all data members which depend on the unfolding results
658}
659
660////////////////////////////////////////////////////////////////////////
661/// Matrix calculations required to propagate systematic errors
663{
664 // data members modified
665 // fEmatUncorrX, fEmatUncorrAx, fDeltaCorrX, fDeltaCorrAx
666 if(!fEmatUncorrX) {
668 }
669 TMatrixDSparse *AM0=nullptr,*AM1=nullptr;
670 if(!fEmatUncorrAx) {
672 if(!AM1) {
674 Int_t *rows_cols=new Int_t[GetNy()];
675 Double_t *data=new Double_t[GetNy()];
676 for(Int_t i=0;i<GetNy();i++) {
677 rows_cols[i]=i;
678 data[i]=1.0;
679 }
682 delete[] data;
683 delete[] rows_cols;
684 AddMSparse(AM1,-1.,one);
687 }
688 }
689 if((!fDeltaSysTau )&&(fDtau>0.0)) {
694 for(Int_t i=0;i<n;i++) {
695 data[i] *= scale;
696 }
697 }
698
700 const TObjString *key;
701
702 // calculate individual systematic errors
703 for(key=(const TObjString *)sysErrIn.Next();key;
704 key=(const TObjString *)sysErrIn.Next()) {
705 const TMatrixDSparse *dsys=
706 (const TMatrixDSparse *)((const TPair *)*sysErrIn)->Value();
707 const TPair *named_emat=(const TPair *)
709 if(!named_emat) {
711 fDeltaCorrX->Add(new TObjString(*key),emat);
712 }
714 if(!named_emat) {
716 if(!AM1) {
718 Int_t *rows_cols=new Int_t[GetNy()];
719 Double_t *data=new Double_t[GetNy()];
720 for(Int_t i=0;i<GetNy();i++) {
721 rows_cols[i]=i;
722 data[i]=1.0;
723 }
726 delete[] data;
727 delete[] rows_cols;
728 AddMSparse(AM1,-1.,one);
731 }
733 fDeltaCorrAx->Add(new TObjString(*key),emat);
734 }
735 }
738}
739
740////////////////////////////////////////////////////////////////////////
741/// Covariance contribution from uncorrelated uncertainties of the
742/// response matrix
743///
744/// \param[inout] ematrix covariance matrix histogram
745/// \param[in] binMap mapping of histogram bins
746/// \param[in] clearEmat if true, ematrix is cleared prior to adding
747/// this covariance matrix contribution
748///
749/// This method propagates the uncertainties of the response matrix
750/// histogram, specified with the constructor, to the unfolding
751/// result. It is assumed that the entries of that histogram are
752/// bin-to-bin uncorrelated. In many cases this corresponds to the
753/// "Monte Carlo statistical uncertainties".
754/// <br/>
755/// The array <b>binMap</b> is explained with the method GetOutput().
756/// The flag <b>clearEmat</b> may be used to add covariance matrices from
757/// several uncertainty sources.
758///
761{
762 // data members modified:
763 // fVYAx, fESparse, fEAtV, fErrorAStat
766}
767
768////////////////////////////////////////////////////////////////////////
769/// propagate uncorrelated systematic errors to a covariance matrix
770///
771/// \param[in] m_0 coefficients for error propagation
772/// \param[in] m_1 coefficients for error propagation
773///
774/// returns the covariance matrix
775
777(const TMatrixDSparse *m_0,const TMatrixDSparse *m_1)
778{
779 // propagate uncorrelated systematic errors to a covariance matrix
780 // m0,m1 : coefficients (matrices) for propagating the errors
781 //
782 // the error matrix is calculated by standard error propagation, where the
783 // derivative of the result vector X wrt the matrix A is given by
784 //
785 // dX_k / dA_ij = M0_kj * Z0_i - M1_ki * Z1_j
786 //
787 // where:
788 // the matrices M0 and M1 are arguments to this function
789 // the vectors Z0, Z1 : GetDXDAZ()
790 //
791 // The matrix A is calculated from a matrix B as
792 //
793 // A_ij = B_ij / sum_k B_kj
794 //
795 // where k runs over additional indices of B, not present in A.
796 // (underflow and overflow bins, used for efficiency corrections)
797 //
798 // define: Norm_j = sum_k B_kj (data member fSumOverY)
799 //
800 // the derivative of A wrt this input matrix B is given by:
801 //
802 // dA_ij / dB_kj = ( delta_ik - A_ij ) * 1/Norm_j
803 //
804 // The covariance matrix Vxx is:
805 //
806 // Vxx_mn = sum_ijlk [ (dX_m / dA_ij) * (dA_ij / dB_kj) * DB_kj
807 // * (dX_n / dA_lj) * (dA_lj / dB_kj) ]
808 //
809 // where DB_kj is the error on B_kj squared
810 // Simplify the sum over k:
811 //
812 // sum_k [ (dA_ij / dB_kj) * DB_kj * (dA_lj / dB_kj) ]
813 // = sum_k [ ( delta_ik - A_ij ) * 1/Norm_j * DB_kj *
814 // * ( delta_lk - A_lj ) * 1/Norm_j ]
815 // = sum_k [ ( delta_ik*delta_lk - delta_ik*A_lj - delta_lk*A_ij
816 // + A_ij * A_lj ) * DB_kj / Norm_j^2 ]
817 //
818 // introduce normalized errors: Rsq_kj = DB_kj / Norm_j^2
819 // after summing over k:
820 // delta_ik*delta_lk*Rsq_kj -> delta_il*Rsq_ij
821 // delta_ik*A_lj*Rsq_kj -> A_lj*Rsq_ij
822 // delta_lk*A_ij*Rsq_kj -> A_ij*Rsq_lj
823 // A_ij*A_lj*Rsq_kj -> A_ij*A_lj*sum_k(Rsq_kj)
824 //
825 // introduce sum of normalized errors squared: SRsq_j = sum_k(Rsq_kj)
826 //
827 // Note: Rsq_ij is stored as fDAinRelSq (excludes extra indices of B)
828 // and SRsq_j is stored as fDAinColRelSq (sum includes all indices of B)
829 //
830 // Vxx_nm = sum_ijl [ (dX_m / dA_ij) * (dX_n / dA_lj)
831 // (delta_il*Rsq_ij - A_lj*Rsq_ij - A_ij*Rsq_lj + A_ij*A_lj *SRsq_j) ]
832 //
833 // Vxx_nm = sum_j [ F_mj * F_nj * SRsq_j
834 // - sum_j [ G_mj * F_nj ]
835 // - sum_j [ F_mj * G_nj ]
836 // + sum_ij [ (dX_m / dA_ij) * (dX_n / dA_lj) * Rsq_ij ]
837 //
838 // where:
839 // F_mj = sum_i [ (dX_m / dA_ij) * A_ij ]
840 // G_mj = sum_i [ (dX_m / dA_ij) * Rsq_ij ]
841 //
842 // In order to avoid explicitly calculating the 3-dimensional tensor
843 // (dX_m/dA_ij) the sums are evaluated further, using
844 // dX_k / dA_ij = M0_kj * Z0_i - M1_ki * Z1_j
845 //
846 // F_mj = M0_mj * (A# Z0)_j - (M1 A)_mj Z1_j
847 // G_mj = M0_mj * (Rsq# Z0)_j - (M1 Rsq)_mj Z1_j
848 //
849 // and
850 //
851 // sum_ij [ (dX_m/dA_ij) * (dX_n/dA_ij) * Rsq_ij ] =
852 // sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
853 // + sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
854 // - sum_i [ M1_mi * H_ni + M1_ni * H_mi]
855 // where:
856 // H_mi = Z0_i * sum_j [ M0_mj * Z1_j * Rsq_ij ]
857 //
858 // collect all contributions:
859 // Vxx_nm = r0 -r1 -r2 +r3 +r4 -r5 -r6
860 // r0 = sum_j [ F_mj * F_nj * SRsq_j ]
861 // r1 = sum_j [ G_mj * F_nj ]
862 // r2 = sum_j [ F_mj * G_nj ]
863 // r3 = sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
864 // r4 = sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
865 // r5 = sum_i [ M1_mi * H_ni ]
866 // r6 = sum_i [ M1_ni * H_mi ]
867
868 //======================================================
869 // calculate contributions containing matrices F and G
870 // r0,r1,r2
871 TMatrixDSparse *r=nullptr;
873 // calculate matrices (M1*A)_mj * Z1_j and (M1*Rsq)_mj * Z1_j
878 // calculate vectors A#*Z0 and Rsq#*Z0
882 //calculate matrix F
883 // F_mj = M0_mj * (A# Z0)_j - (M1 A)_mj Z1_j
886 AddMSparse(F,-1.0,M1A_Z1);
887 //calculate matrix G
888 // G_mj = M0_mj * (Rsq# Z0)_j - (M1 Rsq)_mj Z1_j
891 AddMSparse(G,-1.0,M1Rsq_Z1);
896 // r0 = sum_j [ F_mj * F_nj * SRsq_j ]
898 // r1 = sum_j [ G_mj * F_nj ]
900 // r2 = sum_j [ F_mj * G_nj ]
902 // r = r0-r1-r2
903 AddMSparse(r,-1.0,r1);
904 AddMSparse(r,-1.0,r2);
907 DeleteMatrix(&F);
908 DeleteMatrix(&G);
909 }
910 //======================================================
911 // calculate contribution
912 // sum_ij [ (dX_m/dA_ij) * (dX_n/dA_ij) * Rsq_ij ]
913 // (r3,r4,r5,r6)
914 if(fDAinRelSq) {
915 // (Z0_i)^2
917 const Int_t *Z0sq_rows=Z0sq.GetRowIndexArray();
918 Double_t *Z0sq_data=Z0sq.GetMatrixArray();
919 for(int index=0;index<Z0sq_rows[Z0sq.GetNrows()];index++) {
921 }
922 // Z0sqRsq = sum_i (Z_i)^2 * Rsq_ij
924 // r3 = sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
927
928 // (Z1_j)^2
930 const Int_t *Z1sq_rows=Z1sq.GetRowIndexArray();
931 Double_t *Z1sq_data=Z1sq.GetMatrixArray();
932 for(int index=0;index<Z1sq_rows[Z1sq.GetNrows()];index++) {
934 }
935 // Z1sqRsq = sum_j (Z1_j)^2 * Rsq_ij ]
937 // r4 = sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
940
941 // sum_j [ M0_mj * Z1_j * Rsq_ij ]
944 // H_mi = Z0_i * sum_j [ M0_mj * Z1_j * Rsq_ij ]
946 // r5 = sum_i [ M1_mi * H_ni ]
948 // r6 = sum_i [ H_mi * M1_ni ]
950 DeleteMatrix(&H);
951 // r = r0 -r1 -r2 +r3 +r4 -r5 -r6
952 if(r) {
953 AddMSparse(r,1.0,r3);
955 } else {
956 r=r3;
957 r3=nullptr;
958 }
959 AddMSparse(r,1.0,r4);
960 AddMSparse(r,-1.0,r5);
961 AddMSparse(r,-1.0,r6);
965 }
966 return r;
967}
968
969////////////////////////////////////////////////////////////////////////
970/// propagate correlated systematic shift to an output vector
971///
972/// \param[in] m1 coefficients
973/// \param[in] m2 coeffiicients
974/// \param[in] dsys matrix of correlated shifts from this source
975
977(const TMatrixDSparse *m1,const TMatrixDSparse *m2,const TMatrixDSparse *dsys)
978{
979 // propagate correlated systematic shift to output vector
980 // m1,m2 : coefficients for propagating the errors
981 // dsys : matrix of correlated shifts from this source
982
983 // delta_m =
984 // sum{i,j} {
985 // ((*m1)(m,j) * (*fVYAx)(i) - (*m2)(m,i) * (*fX)(j))*dsys(i,j) }
986 // = sum_j (*m1)(m,j) sum_i dsys(i,j) * (*fVYAx)(i)
987 // - sum_i (*m2)(m,i) sum_j dsys(i,j) * (*fX)(j)
988
995 AddMSparse(delta,-1.0,delta2);
997 return delta;
998}
999
1000////////////////////////////////////////////////////////////////////////
1001/// Specify an uncertainty on tau
1002///
1003/// \param[in] delta_tau new uncertainty on tau
1004///
1005/// The default is to have no uncertyainty on tau.
1007{
1008 // set uncertainty on tau
1011}
1012
1013////////////////////////////////////////////////////////////////////////
1014/// correlated one-sigma shifts correspinding to a given systematic uncertainty
1015///
1016/// \param[out] hist_delta histogram to store shifts
1017/// \param[in] name identifier of the background source
1018/// \param[in] binMap (default=nullptr) remapping of histogram bins
1019///
1020/// returns true if the error source was found.
1021/// <br>
1022/// This method returns the shifts of the unfolding result induced by
1023/// varying the identified systematic source by one sigma.
1024/// <br/>
1025/// the array <b>binMap</b> is explained with the method GetOutput().
1027 const Int_t *binMap)
1028{
1030 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1031 const TMatrixDSparse *delta=nullptr;
1032 if(named_emat) {
1033 delta=(TMatrixDSparse *)named_emat->Value();
1034 }
1036 return delta !=nullptr;
1037}
1038
1039////////////////////////////////////////////////////////////////////////
1040/// correlated one-sigma shifts from background normalisation uncertainty
1041///
1042/// \param[out] hist_delta histogram to store shifts
1043/// \param[in] source identifier of the background source
1044/// \param[in] binMap (default=nullptr) remapping of histogram bins
1045///
1046/// returns true if the background source was found.
1047/// <br>
1048/// This method returns the shifts of the unfolding result induced by
1049/// varying the normalisation of the identified background by one sigma.
1050/// <br/>
1051/// the array <b>binMap</b> is explained with the method GetOutput().
1052
1054(TH1 *hist_delta,const char *source,const Int_t *binMap)
1055{
1058 TMatrixDSparse *dx=nullptr;
1059 if(named_err) {
1060 const TMatrixD *dy=(TMatrixD *)named_err->Value();
1062 }
1064 if(dx!=nullptr) {
1065 DeleteMatrix(&dx);
1066 return kTRUE;
1067 }
1068 return kFALSE;
1069}
1070
1071////////////////////////////////////////////////////////////////////////
1072/// correlated one-sigma shifts from shifting tau
1073///
1074/// \param[out] hist_delta histogram to store shifts
1075/// \param[in] source identifier of the background source
1076/// \param[in] binMap (default=nullptr) remapping of histogram bins
1077///
1078/// returns true if the background source was found.
1079/// <br>
1080/// This method returns the shifts of the unfolding result induced by
1081/// varying the normalisation of the identified background by one sigma.
1082/// <br/>
1083/// the array <b>binMap</b> is explained with the method GetOutput().
1084
1086{
1087 // calculate systematic shift from tau variation
1088 // ematrix: output
1089 // binMap: see method GetEmatrix()
1092 return fDeltaSysTau !=nullptr;
1093}
1094
1095////////////////////////////////////////////////////////////////////////
1096/// covariance contribution from a systematic variation of the
1097/// response matrix
1098///
1099/// \param[inout] ematrix covariance matrix histogram
1100/// \param[in] name identifier of the systematic variation
1101/// \param[in] binMap (default=nullptr) remapping of histogram bins
1102/// \param[in] clearEmat (default=true) if true, clear the histogram
1103/// prior to adding the covariance matrix contribution
1104///
1105/// Returns the covariance matrix contribution from shifting the given
1106/// uncertainty source within one sigma
1107/// <br/>
1108/// the array <b>binMap</b> is explained with the method GetOutput().
1109/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1110/// several uncertainty sources.
1111
1113(TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1114{
1116 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1117 TMatrixDSparse *emat=nullptr;
1118 if(named_emat) {
1119 const TMatrixDSparse *delta=(TMatrixDSparse *)named_emat->Value();
1120 emat=MultiplyMSparseMSparseTranspVector(delta,delta,nullptr);
1121 }
1124}
1125
1126////////////////////////////////////////////////////////////////////////
1127/// covariance contribution from background normalisation uncertainty
1128///
1129/// \param[inout] ematrix output histogram
1130/// \param[in] source identifier of the background source
1131/// \param[in] binMap (default=nullptr) remapping of histogram bins
1132/// \param[in] clearEmat (default=true) if true, clear the histogram
1133/// prior to adding the covariance matrix contribution
1134///
1135/// this method returns the uncertainties on the unfolding result
1136/// arising from the background source <b>source</b> and its normalisation
1137/// uncertainty. See method SubtractBackground() how to set the normalisation uncertainty
1138/// <br/>
1139/// the array <b>binMap</b> is explained with the method GetOutput().
1140/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1141/// several uncertainty sources.
1142
1158
1159////////////////////////////////////////////////////////////////////////
1160/// covariance matrix contribution from error on regularisation
1161/// parameter
1162///
1163/// \param[inout] ematrix output histogram
1164/// \param[in] binMap (default=nullptr) remapping of histogram bins
1165/// \param[in] clearEmat (default=true) if true, clear the histogram
1166///
1167/// this method returns the covariance contributions to the unfolding result
1168/// from the assigned uncertainty on the parameter tau, see method
1169/// SetTauError().
1170/// <br>
1171/// the array <b>binMap</b> is explained with the method GetOutput().
1172/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1173/// several uncertainty sources.
1174
1177{
1178 // calculate error matrix from error in regularisation parameter
1179 // ematrix: output
1180 // binMap: see method GetEmatrix()
1181 // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1183 TMatrixDSparse *emat=nullptr;
1184 if(fDeltaSysTau) {
1186 }
1189}
1190
1191////////////////////////////////////////////////////////////////////////
1192/// covariance matrix contribution from input measurement uncertainties
1193///
1194/// \param[inout] ematrix output histogram
1195/// \param[in] binMap (default=nullptr) remapping of histogram bins
1196/// \param[in] clearEmat (default=true) if true, clear the histogram
1197///
1198/// this method returns the covariance contributions to the unfolding result
1199/// from the uncertainties or covariance of the input
1200/// data. In many cases, these are the "statistical uncertainties".
1201/// <br>
1202/// The array <b>binMap</b> is explained with the method GetOutput().
1203/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1204/// several uncertainty sources.
1205
1211
1212////////////////////////////////////////////////////////////////////////
1213/// covariance contribution from background uncorrelated uncertainty
1214///
1215/// \param[in] ematrix output histogram
1216/// \param[in] source identifier of the background source
1217/// \param[in] binMap (default=nullptr) remapping of histogram bins
1218/// \param[in] clearEmat (default=true) if true, clear the histogram
1219///
1220/// this method returns the covariance contributions to the unfolding result
1221/// arising from the background source <b>source</b> and the uncorrelated
1222/// (background histogram uncertainties). Also see method SubtractBackground()
1223/// <br/>
1224/// the array <b>binMap</b> is explained with the method GetOutput().
1225/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1226/// several uncertainty sources.
1227
1240
1241////////////////////////////////////////////////////////////////////////
1242/// propagate an error matrix on the input vector to the unfolding result
1243///
1244/// \param[in] vyy input error matrix
1245/// \param[inout] ematrix histogram to be updated
1246/// \param[in] binMap mapping of histogram bins
1247/// \param[in] clearEmat if set, clear histogram before adding this
1248/// covariance contribution
1251{
1252 // propagate error matrix vyy to the result
1253 // vyy: error matrix on input data fY
1254 // ematrix: output
1255 // binMap: see method GetEmatrix()
1256 // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1258 TMatrixDSparse *em=nullptr;
1259 if(vyy) {
1263 }
1265 DeleteMatrix(&em);
1266}
1267
1268////////////////////////////////////////////////////////////////////////
1269/// Get total error matrix, summing up all contributions
1270///
1271/// \param[out] ematrix histogram which will be filled
1272/// \param[in] binMap (default=nullptr) remapping of histogram bins
1273///
1274/// the array <b>binMap</b> is explained with the method GetOutput().
1276{
1277 // get total error including statistical error
1278 // ematrix: output
1279 // binMap: see method GetEmatrix()
1280 GetEmatrix(ematrix,binMap); // (stat)+(d)+(e)
1283 const TObject *key;
1284
1285 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1287 ((const TObjString *)key)->GetString(),
1288 binMap,kFALSE); // (b)
1289 }
1291}
1292
1293////////////////////////////////////////////////////////////////////////
1294/// determine total error matrix on the vector Ax
1296{
1298
1299 // errors from input vector and from background subtraction
1301
1302 // uncorrelated systematic error
1303 if(fEmatUncorrAx) {
1305 }
1307 const TObject *key;
1308
1309 // correlated systematic errors
1310 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1311 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1315 }
1316 // error on tau
1317 if(fDeltaSysTau) {
1324 }
1325 return emat_sum;
1326}
1327
1328////////////////////////////////////////////////////////////////////////
1329/// determine total error matrix on the vector x
1331{
1333
1334 // errors from input vector and from background subtraction
1336
1337 // uncorrelated systematic error
1338 if(fEmatUncorrX) {
1340 }
1342 const TObject *key;
1343
1344 // correlated systematic errors
1345 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1346 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1350 }
1351 // error on tau
1352 if(fDeltaSysTau) {
1357 }
1358 return emat_sum;
1359}
1360
1361
1362////////////////////////////////////////////////////////////////////////
1363/// calculate total chi**2 including all systematic errors
1364
1366{
1367
1369
1370 Int_t rank=0;
1374 DeleteMatrix(&v);
1375 Double_t r=0.0;
1376 const Int_t *vdy_rows=vdy->GetRowIndexArray();
1377 const Double_t *vdy_data=vdy->GetMatrixArray();
1378 for(Int_t i=0;i<vdy->GetNrows();i++) {
1379 if(vdy_rows[i+1]>vdy_rows[i]) {
1380 r += vdy_data[vdy_rows[i]] * dy(i,0);
1381 }
1382 }
1383 DeleteMatrix(&vdy);
1385 return r;
1386}
1387
1388////////////////////////////////////////////////////////////////////////
1389/// Get global correlatiocn coefficients, summing up all contributions
1390///
1391/// \param[out] rhoi histogram which will be filled
1392/// \param[in] binMap (default=nullptr) remapping of histogram bins
1393/// \param[out] invEmat (default=nullptr) inverse of error matrix
1394///
1395/// return the global correlation coefficients, including all error
1396/// sources. If <b>invEmat</b> is nonzero, the inverse of the error
1397/// matrix is returned in that histogram
1398/// <br/>
1399/// the array <b>binMap</b> is explained with the method GetOutput().
1401{
1402 // get global correlation coefficients including systematic,statistical,background,tau errors
1403 // rhoi: output histogram
1404 // binMap: for each global bin, indicate in which histogram bin
1405 // to store its content
1406 // invEmat: output histogram for inverse of error matrix
1407 // (pointer may zero if inverse is not requested)
1408 ClearHistogram(rhoi,-1.);
1411
1413}
1414
1415////////////////////////////////////////////////////////////////////////
1416/// scale columns of a matrix by the corresponding rows of a vector
1417///
1418/// \par[inout] m matrix
1419/// \par[in] v vector
1420///
1421/// the entries m<sub>ij</sub> are multiplied by v<sub>j</sub>.
1424{
1425 // scale columns of m by the corresponding rows of v
1426 // input:
1427 // m: pointer to sparse matrix of dimension NxM
1428 // v: pointer to matrix of dimension Mx1
1429 if((m->GetNcols() != v->GetNrows())||(v->GetNcols()!=1)) {
1430 Fatal("ScaleColumnsByVector error",
1431 "matrix cols/vector rows %d!=%d OR vector cols %d !=1\n",
1432 m->GetNcols(),v->GetNrows(),v->GetNcols());
1433 }
1434 const Int_t *rows_m=m->GetRowIndexArray();
1435 const Int_t *cols_m=m->GetColIndexArray();
1436 Double_t *data_m=m->GetMatrixArray();
1437 const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
1438 if(v_sparse) {
1439 const Int_t *rows_v=v_sparse->GetRowIndexArray();
1440 const Double_t *data_v=v_sparse->GetMatrixArray();
1441 for(Int_t i=0;i<m->GetNrows();i++) {
1442 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1445 if(index_v<rows_v[j+1]) {
1447 } else {
1448 data_m[index_m] =0.0;
1449 }
1450 }
1451 }
1452 } else {
1453 for(Int_t i=0;i<m->GetNrows();i++) {
1454 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1456 data_m[index_m] *= (*v)(j,0);
1457 }
1458 }
1459 }
1460}
1461
1462////////////////////////////////////////////////////////////////////////
1463/// map delta to hist_delta, possibly summing up bins
1464///
1465/// \param[out] hist_delta result histogram
1466/// \param[in] delta vector to be mapped to the histogram
1467/// \param[in] binMap mapping of histogram bins
1468///
1469/// grooups of bins of <b>delta</b> are mapped to bins of
1470/// <b>hist_delta</b>. The histogram contents are set to the sum over
1471/// the group of bins. The histogram errors are reset to zero.
1472/// <br/>
1473/// The array <b>binMap</b> is explained with the method GetOutput()
1475(TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap)
1476{
1477 // sum over bins of *delta, as defined in binMap,fXToHist
1478 // hist_delta: histogram to return summed vector
1479 // delta: vector to sum and remap
1480 Int_t nbin=hist_delta->GetNbinsX();
1481 Double_t *c=new Double_t[nbin+2];
1482 for(Int_t i=0;i<nbin+2;i++) {
1483 c[i]=0.0;
1484 }
1485 if(delta) {
1487 const Double_t *delta_data=delta->GetMatrixArray();
1488 const Int_t *delta_rows=delta->GetRowIndexArray();
1489 for(Int_t i=0;i<binMapSize;i++) {
1490 Int_t destBinI=binMap ? binMap[i] : i;
1492 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
1494 if(index<delta_rows[srcBinI+1]) {
1496 }
1497 }
1498 }
1499 }
1500 for(Int_t i=0;i<nbin+2;i++) {
1501 hist_delta->SetBinContent(i,c[i]);
1502 hist_delta->SetBinError(i,0.0);
1503 }
1504 delete[] c;
1505}
1506
1507////////////////////////////////////////////////////////////////////////////////
1508/// Get a new list of all systematic uuncertainty sources.
1509///
1510/// The user is responsible for deleting the list
1512 // get list of names of systematic sources
1513 TSortedList *r=new TSortedList();
1514 TMapIter i(fSysIn);
1515 for(const TObject *key=i.Next();key;key=i.Next()) {
1516 r->Add(((TObjString *)key)->Clone());
1517 }
1518 return r;
1519}
1520
1521////////////////////////////////////////////////////////////////////////////////
1522/// Get a new list of all background sources.
1523///
1524/// The user is responsible for deleting the list
1525/// get list of name of background sources
1526
1528 TSortedList *r=new TSortedList();
1529 TMapIter i(fBgrIn);
1530 for(const TObject *key=i.Next();key;key=i.Next()) {
1531 r->Add(((TObjString *)key)->Clone());
1532 }
1533 return r;
1534}
#define c(i)
Definition RSha256.hxx:101
constexpr Bool_t kFALSE
Definition RtypesCore.h:108
constexpr Bool_t kTRUE
Definition RtypesCore.h:107
ROOT::Detail::TRangeCast< T, true > TRangeDynCast
TRangeDynCast is an adapter class that allows the typed iteration through a TCollection.
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void data
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t index
Option_t Option_t TPoint TPoint const char mode
char name[80]
Definition TGX11.cxx:110
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixT< Double_t > TMatrixD
Definition TMatrixDfwd.h:23
Int_t GetSize() const
Definition TArray.h:47
virtual Int_t GetEntries() const
TH1 is the base class of all histogram classes in ROOT.
Definition TH1.h:109
Service class for 2-D histogram classes.
Definition TH2.h:30
Iterator of map.
Definition TMap.h:144
TObject * Next() override
Returns the next key from a map.
Definition TMap.cxx:540
TMap implements an associative array of (key,value) pairs using a THashTable for efficient retrieval ...
Definition TMap.h:40
void Add(TObject *obj) override
This function may not be used (but we need to provide it since it is a pure virtual in TCollection).
Definition TMap.cxx:53
virtual void SetOwnerKeyValue(Bool_t ownkeys=kTRUE, Bool_t ownvals=kTRUE)
Set ownership for keys and values.
Definition TMap.cxx:351
TObject * FindObject(const char *keyname) const override
Check if a (key,value) pair exists with keyname as name of the key.
Definition TMap.cxx:214
void Clear(Option_t *option="") override
Remove all (key,value) pairs from the map.
Definition TMap.cxx:96
Int_t GetNrows() const
const Int_t * GetRowIndexArray() const override
const Int_t * GetColIndexArray() const override
const Element * GetMatrixArray() const override
Collectable string class.
Definition TObjString.h:28
const TString & GetString() const
Definition TObjString.h:46
Mother of all ROOT objects.
Definition TObject.h:41
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition TObject.cxx:1057
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition TObject.cxx:1071
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition TObject.cxx:1099
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition TObject.cxx:1045
Class used by TMap to store (key,value) pairs.
Definition TMap.h:102
A sorted doubly linked list.
Definition TSortedList.h:28
Basic string class.
Definition TString.h:138
TMatrixD * fAoutside
Input: underflow/overflow bins.
Definition TUnfoldSys.h:68
TMatrixDSparse * fDAinRelSq
Input: normalized errors from input matrix.
Definition TUnfoldSys.h:64
TMatrixDSparse * GetSummedErrorMatrixXX(void)
determine total error matrix on the vector x
void GetEmatrixSysTau(TH2 *ematrix, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
covariance matrix contribution from error on regularisation parameter
Double_t GetChi2Sys(void)
calculate total chi**2 including all systematic errors
void GetEmatrixTotal(TH2 *ematrix, const Int_t *binMap=nullptr)
Get total error matrix, summing up all contributions.
void VectorMapToHist(TH1 *hist_delta, const TMatrixDSparse *delta, const Int_t *binMap)
map delta to hist_delta, possibly summing up bins
void ScaleColumnsByVector(TMatrixDSparse *m, const TMatrixTBase< Double_t > *v) const
scale columns of a matrix by the corresponding rows of a vector
void GetEmatrixSysBackgroundScale(TH2 *ematrix, const char *source, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
covariance contribution from background normalisation uncertainty
TMap * fDeltaCorrAx
Result: syst.shift from fSysIn on fAx.
Definition TUnfoldSys.h:90
TMatrixD * fYData
Input: fY prior to bgr subtraction.
Definition TUnfoldSys.h:80
void GetRhoItotal(TH1 *rhoi, const Int_t *binMap=nullptr, TH2 *invEmat=nullptr)
Get global correlatiocn coefficients, summing up all contributions.
void GetEmatrixFromVyy(const TMatrixDSparse *vyy, TH2 *ematrix, const Int_t *binMap, Bool_t clearEmat)
propagate an error matrix on the input vector to the unfolding result
void GetBackground(TH1 *bgr, const char *bgrSource=nullptr, const Int_t *binMap=nullptr, Int_t includeError=3, Bool_t clearHist=kTRUE) const
get background into a histogram
void InitTUnfoldSys(void)
TMatrixDSparse * fVyyData
Input: error on fY prior to bgr subtraction.
Definition TUnfoldSys.h:82
TMatrixDSparse * fEmatUncorrAx
Result: syst.error from fDA2 on fAx.
Definition TUnfoldSys.h:86
Bool_t GetDeltaSysBackgroundScale(TH1 *delta, const char *source, const Int_t *binMap=nullptr)
correlated one-sigma shifts from background normalisation uncertainty
void DoBackgroundSubtraction(void)
perform background subtraction
TMatrixDSparse * GetSummedErrorMatrixYY(void)
determine total error matrix on the vector Ax
Double_t fDtau
Input: error on tau.
Definition TUnfoldSys.h:78
ESysErrMode
type of matrix specified with AddSysError()
Definition TUnfoldSys.h:106
@ kSysErrModeRelative
matrix gives the relative shifts
Definition TUnfoldSys.h:112
@ kSysErrModeMatrix
matrix is an alternative to the default matrix, the errors are the difference to the original matrix
Definition TUnfoldSys.h:108
@ kSysErrModeShift
matrix gives the absolute shifts
Definition TUnfoldSys.h:110
virtual TMatrixDSparse * PrepareUncorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2)
propagate uncorrelated systematic errors to a covariance matrix
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=nullptr)
correlated one-sigma shifts from shifting tau
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=nullptr)
correlated one-sigma shifts correspinding to a given systematic uncertainty
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
covariance matrix contribution from input measurement uncertainties
TMatrixDSparse * fEmatUncorrX
Result: syst.error from fDA2 on fX.
Definition TUnfoldSys.h:84
~TUnfoldSys(void) override
void SubtractBackground(const TH1 *hist_bgr, const char *name, Double_t scale=1.0, Double_t scale_error=0.0)
Specify a source of background.
void ClearResults(void) override
Clear all data members which depend on the unfolding results.
TMap * fBgrIn
Input: size of background sources.
Definition TUnfoldSys.h:72
TMap * fDeltaCorrX
Result: syst.shift from fSysIn on fX.
Definition TUnfoldSys.h:88
TSortedList * GetSysSources(void) const
Get a new list of all systematic uuncertainty sources.
TMatrixD * fDAinColRelSq
Input: normalized column err.sq. (inp.matr.)
Definition TUnfoldSys.h:66
TMatrixDSparse * fDeltaSysTau
Result: systematic shift from tau.
Definition TUnfoldSys.h:92
void AddSysError(const TH2 *sysError, const char *name, EHistMap histmap, ESysErrMode mode)
Specify a correlated systematic uncertainty.
void SetTauError(Double_t delta_tau)
Specify an uncertainty on tau.
TMap * fBgrErrUncorrInSq
Input: uncorr error squared from bgr sources.
Definition TUnfoldSys.h:74
TUnfoldSys(void)
only for use by root streamer or derived classes
virtual void PrepareSysError(void)
Matrix calculations required to propagate systematic errors.
Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=nullptr, const TH2 *hist_vyy_inv=nullptr) override
Define input data for subsequent calls to DoUnfold(tau)
void GetEmatrixSysSource(TH2 *ematrix, const char *source, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
covariance contribution from a systematic variation of the response matrix
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
covariance contribution from background uncorrelated uncertainty
void GetEmatrixSysUncorr(TH2 *ematrix, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
Covariance contribution from uncorrelated uncertainties of the response matrix.
TMap * fBgrErrScaleIn
Input: background sources correlated error.
Definition TUnfoldSys.h:76
TSortedList * GetBgrSources(void) const
Get a new list of all background sources.
virtual TMatrixDSparse * PrepareCorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixDSparse *dsys)
propagate correlated systematic shift to an output vector
TMap * fSysIn
Input: correlated errors.
Definition TUnfoldSys.h:70
An algorithm to unfold distributions from detector to truth level.
Definition TUnfold.h:107
TArrayI fHistToX
mapping of histogram bins to matrix indices
Definition TUnfold.h:170
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse *a, const TMatrixD *b) const
multiply sparse matrix and a non-sparse matrix
Definition TUnfold.cxx:760
TMatrixDSparse * MultiplyMSparseTranspMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
multiply a transposed Sparse matrix with another Sparse matrix
Definition TUnfold.cxx:678
TMatrixDSparse * MultiplyMSparseMSparseTranspVector(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixTBase< Double_t > *v) const
calculate a sparse matrix product M1*V*M2T where the diagonal matrix V is given by a vector
Definition TUnfold.cxx:820
TMatrixDSparse * CreateSparseMatrix(Int_t nrow, Int_t ncol, Int_t nele, Int_t *row, Int_t *col, Double_t *data) const
create a sparse matrix, given the nonzero elements
Definition TUnfold.cxx:579
const TMatrixDSparse * GetDXDAM(int i) const
matrix contributions of the derivative dx/dA
Definition TUnfold.h:252
Int_t GetNy(void) const
returns the number of measurement bins
Definition TUnfold.h:238
const TMatrixDSparse * GetDXDtauSquared(void) const
vector of derivative dx/dtauSquared, using internal bin counting
Definition TUnfold.h:265
static void DeleteMatrix(TMatrixD **m)
delete matrix and invalidate pointer
Definition TUnfold.cxx:188
void ClearHistogram(TH1 *h, Double_t x=0.) const
Initialize bin contents and bin errors for a given histogram.
Definition TUnfold.cxx:3680
Int_t GetNx(void) const
returns internal number of output (truth) matrix rows
Definition TUnfold.h:230
const TMatrixDSparse * GetDXDAZ(int i) const
vector contributions of the derivative dx/dA
Definition TUnfold.h:254
EConstraint
type of extra constraint
Definition TUnfold.h:113
TMatrixDSparse * fVyy
covariance matrix Vyy corresponding to y
Definition TUnfold.h:160
const TMatrixDSparse * GetVyyInv(void) const
inverse of covariance matrix of the data y
Definition TUnfold.h:260
TArrayD fSumOverY
truth vector calculated from the non-normalized response matrix
Definition TUnfold.h:172
ERegMode
choice of regularisation scheme
Definition TUnfold.h:123
void ErrorMatrixToHist(TH2 *ematrix, const TMatrixDSparse *emat, const Int_t *binMap, Bool_t doClear) const
add up an error matrix, also respecting the bin mapping
Definition TUnfold.cxx:3379
TArrayI fXToHist
mapping of matrix indices to histogram bins
Definition TUnfold.h:168
const TMatrixDSparse * GetAx(void) const
vector of folded-back result
Definition TUnfold.h:248
TMatrixD * fY
input (measured) data y
Definition TUnfold.h:158
TMatrixDSparse * InvertMSparseSymmPos(const TMatrixDSparse *A, Int_t *rank) const
get the inverse or pseudo-inverse of a positive, sparse matrix
Definition TUnfold.cxx:993
Double_t fTauSquared
regularisation parameter tau squared
Definition TUnfold.h:166
virtual void ClearResults(void)
reset all results
Definition TUnfold.cxx:208
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=nullptr) const
get output covariance matrix, possibly cumulated over several bins
Definition TUnfold.cxx:3446
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
multiply two sparse matrices
Definition TUnfold.cxx:603
EHistMap
arrangement of axes for the response matrix (TH2 histogram)
Definition TUnfold.h:143
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition TUnfold.h:146
void AddMSparse(TMatrixDSparse *dest, Double_t f, const TMatrixDSparse *src) const
add a sparse matrix, scaled by a factor, to another scaled matrix
Definition TUnfold.cxx:915
const TMatrixDSparse * GetVxx(void) const
covariance matrix of the result
Definition TUnfold.h:244
Double_t GetRhoIFromMatrix(TH1 *rhoi, const TMatrixDSparse *eOrig, const Int_t *binMap, TH2 *invEmat) const
Definition TUnfold.cxx:3553
const TMatrixDSparse * GetDXDY(void) const
matrix of derivatives dx/dy
Definition TUnfold.h:250
TMatrixDSparse * fA
response matrix A
Definition TUnfold.h:154
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=nullptr, const TH2 *hist_vyy_inv=nullptr)
Define input data for subsequent calls to DoUnfold(tau)
Definition TUnfold.cxx:2274
const Int_t n
Definition legend1.C:16
#define G(x, y, z)
#define H(x, y, z)
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition TMath.h:673
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition TMath.h:732
TMarker m
Definition textangle.C:8
static uint64_t sum(uint64_t i)
Definition Factory.cxx:2339