Logo ROOT  
Reference Guide
TUnfoldSys.cxx
Go to the documentation of this file.
1// @(#)root/unfold:$Id$
2// Author: Stefan Schmitt DESY, 23/01/09
3
4/** \class TUnfoldSys
5\ingroup Unfold
6An algorithm to unfold distributions from detector to truth level,
7with background subtraction and propagation of systematic uncertainties
8
9TUnfoldSys is used to decompose a measurement y into several sources x,
10given the measurement uncertainties, background b and a matrix of migrations A.
11The method can be applied to a large number of problems,
12where the measured distribution y is a linear superposition
13of several Monte Carlo shapes. Beyond such a simple template fit,
14TUnfoldSys has an adjustable regularisation term and also supports an
15optional constraint on the total number of events.
16Background sources can be specified, with a normalisation constant and
17normalisation uncertainty. In addition, variants of the response
18matrix may be specified, these are taken to determine systematic
19uncertainties.
20
21<b>For most applications, it is better to use the derived class
22TUnfoldDensity instead of TUnfoldSys. TUnfoldDensity adds
23features to TUnfoldSys, related to possible complex multidimensional
24arrangements of bins. For innocent
25users, the most notable improvement of TUnfoldDensity over TUnfoldSys are
26the getter functions. For TUnfoldSys, histograms have to be booked by the
27user and the getter functions fill the histogram bins. TUnfoldDensity
28simply returns a new, already filled histogram.</b>
29
30If you use this software, please consider the following citation
31
32<b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
33
34Detailed documentation and updates are available on
35http://www.desy.de/~sschmitt
36
37Brief recipy to use TUnfoldSys:
38
39 - a matrix (truth,reconstructed) is given as a two-dimensional histogram
40 as argument to the constructor of TUnfold
41 - a vector of measurements is given as one-dimensional histogram using
42 the SetInput() method
43 - repeated calls to SubtractBackground() to specify background sources
44 - repeated calls to AddSysError() to specify systematic uncertainties
45 - The unfolding is performed
46 - either once with a fixed parameter tau, method DoUnfold(tau)
47 - or multiple times in a scan to determine the best chouce of tau,
48 method ScanLCurve()
49 - Unfolding results are retrieved using various GetXXX() methods
50
51
52Description of (systematic) uncertainties available in
53TUnfoldSys. There are covariance matrix contributions and there are
54systematic shifts. Systematic shifts correspond to the variation of a
55(buicance) parameter, for example a background normalisation or a
56one-sigma variation of a correlated systematic error.
57
58| | Set by | Access covariance matrix | Access vector of shifts | Description |
59|-------------------------|------------------------|---------------------------------|------------------------------|-------------|
60| (a) | TUnfoldSys constructor | GetEmatrixSysUncorr() | n.a. | uncorrelated errors on the input matrix histA, taken as the errors provided with the histogram. These are typically statistical errors from finite Monte Carlo samples. |
61| (b) | AddSysError() | GetEmatrixSysSource() | GetDeltaSysSource() | correlated shifts of the input matrix histA. These shifts are taken as one-sigma effects when switchig on a given error soure. Several such error sources may be defined |
62| (c) | SetTauError() | GetEmatrixSysTau() | GetDeltaSysTau() | A systematic error on the regularisation parameter tau |
63| (d) | SubtractBackground() | GetEmatrixSysBackgroundUncorr() | n.a. | uncorrelated errors on background sources, originating from the errors provided with the background histograms |
64| (e) | SubtractBackground() | GetEmatrixSysBackgroundScale() | GetDeltaSysBackgroundScale() | scale errors on background sources |
65| (i) | SetInput() | GetEmatrixInput() | n.a. | statistical uncertainty of the input (the measurement) |
66| (i)+(d)+(e) | see above | GetEmatrix() | n.a. | Partial sun of uncertainties: all sources which are propagated to the covariance before unfolding |
67| (i)+(a)+(b)+(c)+(d)+(e) | see above | GetEmatrixTotal() | n.a. | All known error sources summed up |
68
69
70Note: (a), (b), (c) are propagated to the result AFTER unfolding,
71whereas the background errors (d) and (e) are added to the data errors
72BEFORE unfolding. For this reason the errors of type (d) and (e) are
73INCLUDED in the standard error matrix and other methods provided by
74the base class TUnfold, whereas errors of type (a), (b), (c) are NOT
75INCLUDED in the methods provided by the base class TUnfold.
76
77
78--------------------------------------------------------------------------------
79<b>Version 17.6, with updated doxygen comments</b>
80
81#### History:
82 - Version 17.5, in parallel to changes in TUnfold
83 - Version 17.4, in parallel to changes in TUnfoldBinning
84 - Version 17.3, in parallel to changes in TUnfoldBinning
85 - Version 17.2, add methods to find back systematic and background sources
86 - Version 17.1, bug fix with background uncertainty
87 - Version 17.0, possibility to specify an error matrix with SetInput
88 - Version 16.1, parallel to changes in TUnfold
89 - Version 16.0, parallel to changes in TUnfold
90 - Version 15, fix bugs with uncorr. uncertainties, add backgnd subtraction
91 - Version 14, remove some print-out, do not add unused sys.errors
92 - Version 13, support for systematic errors This file is part of TUnfold.
93
94 TUnfold is free software: you can redistribute it and/or modify
95 it under the terms of the GNU General Public License as published by
96 the Free Software Foundation, either version 3 of the License, or
97 (at your option) any later version.
98
99 TUnfold is distributed in the hope that it will be useful,
100 but WITHOUT ANY WARRANTY; without even the implied warranty of
101 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
102 GNU General Public License for more details.
103
104 You should have received a copy of the GNU General Public License
105 along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
106*/
107
108#include <iostream>
109#include <TMap.h>
110#include <TMath.h>
111#include <TObjString.h>
112#include <TSortedList.h>
113#include <cmath>
114
115#include "TUnfoldSys.h"
116
118
120{
121 // delete all data members
124 delete fBgrIn;
125 delete fBgrErrUncorrInSq;
126 delete fBgrErrScaleIn;
127 delete fSysIn;
128 ClearResults();
129 delete fDeltaCorrX;
130 delete fDeltaCorrAx;
133}
134
135////////////////////////////////////////////////////////////////////////////////
136/// Only for use by root streamer or derived classes.
137
139{
140 // set all pointers to zero
142}
143
144////////////////////////////////////////////////////////////////////////////////
145/// Set up response matrix A, uncorrelated uncertainties of A and
146/// regularisation scheme.
147///
148/// \param[in] hist_A matrix that describes the migrations
149/// \param[in] histmap mapping of the histogram axes to the unfolding output
150/// \param[in] regmode (default=kRegModeSize) global regularisation mode
151/// \param[in] constraint (default=kEConstraintArea) type of constraint
152///
153/// For further details, consult the constructir of the class TUnfold.
154/// The uncertainties of hist_A are taken to be uncorrelated and aper
155/// propagated to the unfolding result, method GetEmatrixSysUncorr().
156
158(const TH2 *hist_A, EHistMap histmap, ERegMode regmode,EConstraint constraint)
159 : TUnfold(hist_A,histmap,regmode,constraint)
160{
161 // data members initialized to something different from zero:
162 // fDA2, fDAcol
163
164 // initialize TUnfoldSys
166
167 // save underflow and overflow bins
168 fAoutside = new TMatrixD(GetNx(),2);
169 // save the normalized errors on hist_A
170 // to the matrices fDAinRelSq and fDAinColRelSq
171 fDAinColRelSq = new TMatrixD(GetNx(),1);
172
173 Int_t nmax=GetNx()*GetNy();
174 Int_t *rowDAinRelSq = new Int_t[nmax];
175 Int_t *colDAinRelSq = new Int_t[nmax];
176 Double_t *dataDAinRelSq = new Double_t[nmax];
177
178 Int_t da_nonzero=0;
179 for (Int_t ix = 0; ix < GetNx(); ix++) {
180 Int_t ibinx = fXToHist[ix];
181 Double_t sum_sq= fSumOverY[ix]*fSumOverY[ix];
182 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
183 Double_t dz;
184 if (histmap == kHistMapOutputHoriz) {
185 dz = hist_A->GetBinError(ibinx, ibiny);
186 } else {
187 dz = hist_A->GetBinError(ibiny, ibinx);
188 }
189 Double_t normerr_sq=dz*dz/sum_sq;
190 // quadratic sum of all errors from all bins,
191 // including under/overflow bins
192 (*fDAinColRelSq)(ix,0) += normerr_sq;
193
194 if(ibiny==0) {
195 // underflow bin
196 if (histmap == kHistMapOutputHoriz) {
197 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibinx, ibiny);
198 } else {
199 (*fAoutside)(ix,0)=hist_A->GetBinContent(ibiny, ibinx);
200 }
201 } else if(ibiny==GetNy()+1) {
202 // overflow bins
203 if (histmap == kHistMapOutputHoriz) {
204 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibinx, ibiny);
205 } else {
206 (*fAoutside)(ix,1)=hist_A->GetBinContent(ibiny, ibinx);
207 }
208 } else {
209 // error on this bin
210 rowDAinRelSq[da_nonzero]=ibiny-1;
211 colDAinRelSq[da_nonzero] = ix;
212 dataDAinRelSq[da_nonzero] = normerr_sq;
213 if(dataDAinRelSq[da_nonzero]>0.0) da_nonzero++;
214 }
215 }
216 }
217 if(da_nonzero) {
218 fDAinRelSq = CreateSparseMatrix(GetNy(),GetNx(),da_nonzero,
219 rowDAinRelSq,colDAinRelSq,dataDAinRelSq);
220 } else {
222 }
223 delete[] rowDAinRelSq;
224 delete[] colDAinRelSq;
225 delete[] dataDAinRelSq;
226}
227
228////////////////////////////////////////////////////////////////////////////////
229/// Specify a correlated systematic uncertainty.
230///
231/// \param[in] sysError alternative matrix or matrix of absolute/relative shifts
232/// \param[in] name identifier of the error source
233/// \param[in] histmap mapping of the histogram axes
234/// \param[in] mode format of the error source
235///
236/// <b>sysError</b> corresponds to a one-sigma variation. If
237/// may be given in various forms, specified by <b>mode</b>
238///
239/// - <b>mode=kSysErrModeMatrix</b> the histogram <b>sysError</b>
240/// corresponds to an alternative response matrix.
241/// - <b>mode=kSysErrModeShift</b> the content of the histogram <b>sysError</b> are the absolute shifts of the response matrix
242/// - <b>mode=kSysErrModeRelative</b> the content of the histogram <b>sysError</b>
243/// specifies the relative uncertainties
244///
245/// Internally, all three cases are transformed to the case <b>mode=kSysErrModeMatrix</b>.
246
248(const TH2 *sysError,const char *name,EHistMap histmap,ESysErrMode mode)
249{
250
251 if(fSysIn->FindObject(name)) {
252 Error("AddSysError","Source %s given twice, ignoring 2nd call.\n",name);
253 } else {
254 // a copy of fA is made. It can be accessed inside the loop
255 // without having to take care that the sparse structure of *fA
256 // otherwise, *fA may be accidentally destroyed by asking
257 // for an element which is zero.
258 TMatrixD aCopy(*fA);
259
260 Int_t nmax= GetNx()*GetNy();
261 Double_t *data=new Double_t[nmax];
262 Int_t *cols=new Int_t[nmax];
263 Int_t *rows=new Int_t[nmax];
264 nmax=0;
265 for (Int_t ix = 0; ix < GetNx(); ix++) {
266 Int_t ibinx = fXToHist[ix];
267 Double_t sum=0.0;
268 for(Int_t loop=0;loop<2;loop++) {
269 for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
270 Double_t z;
271 // get bin content, depending on histmap
272 if (histmap == kHistMapOutputHoriz) {
273 z = sysError->GetBinContent(ibinx, ibiny);
274 } else {
275 z = sysError->GetBinContent(ibiny, ibinx);
276 }
277 // correct to absolute numbers
279 Double_t z0;
280 if((ibiny>0)&&(ibiny<=GetNy())) {
281 z0=aCopy(ibiny-1,ix)*fSumOverY[ix];
282 } else if(ibiny==0) {
283 z0=(*fAoutside)(ix,0);
284 } else {
285 z0=(*fAoutside)(ix,1);
286 }
288 z += z0;
289 } else if(mode==kSysErrModeRelative) {
290 z = z0*(1.+z);
291 }
292 }
293 if(loop==0) {
294 // sum up all entries, including overflow bins
295 sum += z;
296 } else {
297 if((ibiny>0)&&(ibiny<=GetNy())) {
298 // save normalized matrix of shifts,
299 // excluding overflow bins
300 rows[nmax]=ibiny-1;
301 cols[nmax]=ix;
302 if(sum>0.0) {
303 data[nmax]=z/sum-aCopy(ibiny-1,ix);
304 } else {
305 data[nmax]=0.0;
306 }
307 if(data[nmax] != 0.0) nmax++;
308 }
309 }
310 }
311 }
312 }
313 if(nmax==0) {
314 Error("AddSysError",
315 "source %s has no influence and has not been added.\n",name);
316 } else {
318 nmax,rows,cols,data);
319 fSysIn->Add(new TObjString(name),dsys);
320 }
321 delete[] data;
322 delete[] rows;
323 delete[] cols;
324 }
325}
326
327////////////////////////////////////////////////////////////////////////////////
328/// Perform background subtraction.
329///
330/// This prepares the data members for the base class TUnfold, such
331/// that the background is properly taken into account.
332
334{
335 // fY = fYData - fBgrIn
336 // fVyy = fVyyData + fBgrErrUncorr^2 + fBgrErrCorr * fBgrErrCorr#
337 // fVyyinv = fVyy^(-1)
338
339 if(fYData) {
342 if(fBgrIn->GetEntries()<=0) {
343 // simple copy
344 fY=new TMatrixD(*fYData);
346 } else {
347 if(GetVyyInv()) {
348 Warning("DoBackgroundSubtraction",
349 "inverse error matrix from user input,"
350 " not corrected for background");
351 }
352 // copy of the data
353 fY=new TMatrixD(*fYData);
354 // subtract background from fY
355 const TObject *key;
356 {
357 TMapIter bgrPtr(fBgrIn);
358 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
359 const TMatrixD *bgr=(const TMatrixD *)((const TPair *)*bgrPtr)->Value();
360 for(Int_t i=0;i<GetNy();i++) {
361 (*fY)(i,0) -= (*bgr)(i,0);
362 }
363 }
364 }
365 // copy original error matrix
366 TMatrixD vyy(*fVyyData);
367 // determine used bins
368 Int_t ny=fVyyData->GetNrows();
369 const Int_t *vyydata_rows=fVyyData->GetRowIndexArray();
370 const Int_t *vyydata_cols=fVyyData->GetColIndexArray();
371 const Double_t *vyydata_data=fVyyData->GetMatrixArray();
372 Int_t *usedBin=new Int_t[ny];
373 for(Int_t i=0;i<ny;i++) {
374 usedBin[i]=0;
375 }
376 for(Int_t i=0;i<ny;i++) {
377 for(Int_t k=vyydata_rows[i];k<vyydata_rows[i+1];k++) {
378 if(vyydata_data[k]>0.0) {
379 usedBin[i]++;
380 usedBin[vyydata_cols[k]]++;
381 }
382 }
383 }
384 // add uncorrelated background errors
385 {
386 TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
387 for(key=bgrErrUncorrSqPtr.Next();key;
388 key=bgrErrUncorrSqPtr.Next()) {
389 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)((const TPair *)*bgrErrUncorrSqPtr)->Value();
390 for(Int_t yi=0;yi<ny;yi++) {
391 if(!usedBin[yi]) continue;
392 vyy(yi,yi) +=(*bgrerruncorrSquared)(yi,0);
393 }
394 }
395 }
396 // add correlated background errors
397 {
398 TMapIter bgrErrScalePtr(fBgrErrScaleIn);
399 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
400 const TMatrixD *bgrerrscale=(const TMatrixD *)((const TPair *)*bgrErrScalePtr)->Value();
401 for(Int_t yi=0;yi<ny;yi++) {
402 if(!usedBin[yi]) continue;
403 for(Int_t yj=0;yj<ny;yj++) {
404 if(!usedBin[yj]) continue;
405 vyy(yi,yj) +=(*bgrerrscale)(yi,0)* (*bgrerrscale)(yj,0);
406 }
407 }
408 }
409 }
410 delete[] usedBin;
411 usedBin=0;
412
413 // convert to sparse matrix
414 fVyy=new TMatrixDSparse(vyy);
415
416 }
417 } else {
418 Fatal("DoBackgroundSubtraction","No input vector defined");
419 }
420}
421
422////////////////////////////////////////////////////////////////////////////////
423/// Define the input data for subsequent calls to DoUnfold(Double_t).
424///
425/// input: input distribution with errors
426/// - scaleBias: scale factor applied to the bias
427/// - oneOverZeroError: for bins with zero error, this number defines 1/error.
428///
429/// Return value: number of bins with bad error
430/// +10000*number of unconstrained output bins
431///
432/// Note: return values>=10000 are fatal errors,
433/// for the given input, the unfolding can not be done!
434///
435/// Calls the SetInput method of the base class, then renames the input
436/// vectors fY and fVyy, then performs the background subtraction
437///
438/// Data members modified:
439/// fYData,fY,fVyyData,fVyy,fVyyinvData,fVyyinv
440///
441/// and those modified by TUnfold::SetInput()
442/// and those modified by DoBackgroundSubtraction()
443
444Int_t TUnfoldSys::SetInput(const TH1 *hist_y,Double_t scaleBias,
445 Double_t oneOverZeroError,const TH2 *hist_vyy,
446 const TH2 *hist_vyy_inv)
447{
448
449 Int_t r=TUnfold::SetInput(hist_y,scaleBias,oneOverZeroError,hist_vyy,
450 hist_vyy_inv);
451 fYData=fY;
452 fY=0;
454 fVyy=0;
456
457 return r;
458}
459
460////////////////////////////////////////////////////////////////////////////////
461/// Specify a source of background.
462///
463/// \param[in] bgr background distribution with uncorrelated errors
464/// \param[in] name identifier for this background source
465/// \param[in] scale normalisation factor applied to the background
466/// \param[in] scale_error normalisation uncertainty
467///
468/// The contribution <b>scale</b>*<b>bgr</b> is subtracted from the
469/// measurement prior to unfolding. The following contributions are
470/// added to the input covarianc ematrix
471///
472/// - using the uncorrelated histogram errors <b>dbgr</b>, the contribution
473/// (<b>scale</b>*<b>dbgr<sub>i</sub></b>)<sup>2</sup> is added to the
474/// diagonals of the covariance
475/// - using the histogram contents, the background normalisation uncertainty contribution
476/// <b>dscale</b>*<b>bgr<sub>i</sub></b> <b>dscale</b>*<b>bgr<sub>j</sub></b>
477/// is added to the covariance matrix
478///
479/// Data members modified:
480/// fBgrIn,fBgrErrUncorrInSq,fBgrErrScaleIn and those modified by DoBackgroundSubtraction()
481
483(const TH1 *bgr,const char *name,Double_t scale,Double_t scale_error)
484{
485
486 // save background source
487 if(fBgrIn->FindObject(name)) {
488 Error("SubtractBackground","Source %s given twice, ignoring 2nd call.\n",
489 name);
490 } else {
491 TMatrixD *bgrScaled=new TMatrixD(GetNy(),1);
492 TMatrixD *bgrErrUncSq=new TMatrixD(GetNy(),1);
493 TMatrixD *bgrErrCorr=new TMatrixD(GetNy(),1);
494 for(Int_t row=0;row<GetNy();row++) {
495 (*bgrScaled)(row,0) = scale*bgr->GetBinContent(row+1);
496 (*bgrErrUncSq)(row,0) =
497 TMath::Power(scale*bgr->GetBinError(row+1),2.);
498 (*bgrErrCorr)(row,0) = scale_error*bgr->GetBinContent(row+1);
499 }
500 fBgrIn->Add(new TObjString(name),bgrScaled);
501 fBgrErrUncorrInSq->Add(new TObjString(name),bgrErrUncSq);
502 fBgrErrScaleIn->Add(new TObjString(name),bgrErrCorr);
503 if(fYData) {
505 } else {
506 Info("SubtractBackground",
507 "Background subtraction prior to setting input data");
508 }
509 }
510}
511
512////////////////////////////////////////////////////////////////////////////////
513/// Get background into a histogram.
514///
515/// \param[inout] bgrHist target histogram, content and errors will be altered
516/// \param[in] bgrSource (default=0) name of backgrond source or zero
517/// to add all sources of background
518/// \param[in] binMap (default=0) remap histogram bins
519/// \param[in] includeError (default=3) include uncorrelated(1),
520/// correlated (2) or both (3) sources of uncertainty in the
521/// histogram errors
522/// \param[in] clearHist (default=true) reset histogram before adding
523/// up the specified background sources
524///
525/// the array <b>binMap</b> is explained with the method GetOutput().
526/// The flag <b>clearHist</b> may be used to add background from
527/// several sources in successive calls to GetBackground().
528
530(TH1 *bgrHist,const char *bgrSource,const Int_t *binMap,
531 Int_t includeError,Bool_t clearHist) const
532{
533 if(clearHist) {
534 ClearHistogram(bgrHist);
535 }
536 // add all background sources
537 const TObject *key;
538 {
539 TMapIter bgrPtr(fBgrIn);
540 for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
541 TString bgrName=((const TObjString *)key)->GetString();
542 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
543 const TMatrixD *bgr=(const TMatrixD *)((const TPair *)*bgrPtr)->Value();
544 for(Int_t i=0;i<GetNy();i++) {
545 Int_t destBin=binMap[i];
546 bgrHist->SetBinContent(destBin,bgrHist->GetBinContent(destBin)+
547 (*bgr)(i,0));
548 }
549 }
550 }
551 // add uncorrelated background errors
552 if(includeError &1) {
553 TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
554 for(key=bgrErrUncorrSqPtr.Next();key;key=bgrErrUncorrSqPtr.Next()) {
555 TString bgrName=((const TObjString *)key)->GetString();
556 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
557 const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)
558 ((const TPair *)*bgrErrUncorrSqPtr)->Value();
559 for(Int_t i=0;i<GetNy();i++) {
560 Int_t destBin=binMap[i];
561 bgrHist->SetBinError
562 (destBin,TMath::Sqrt
563 ((*bgrerruncorrSquared)(i,0)+
564 TMath::Power(bgrHist->GetBinError(destBin),2.)));
565 }
566 }
567 }
568 if(includeError & 2) {
569 TMapIter bgrErrScalePtr(fBgrErrScaleIn);
570 for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
571 TString bgrName=((const TObjString *)key)->GetString();
572 if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
573 const TMatrixD *bgrerrscale=(TMatrixD const *)((const TPair *)*bgrErrScalePtr)->Value();
574 for(Int_t i=0;i<GetNy();i++) {
575 Int_t destBin=binMap[i];
576 bgrHist->SetBinError(destBin,hypot((*bgrerrscale)(i,0),
577 bgrHist->GetBinError(destBin)));
578 }
579 }
580 }
581}
582
583////////////////////////////////////////////////////////////////////////////////
584/// Initialize pointers and TMaps.
585
587{
588 // input
589 fDAinRelSq = 0;
590 fDAinColRelSq = 0;
591 fAoutside = 0;
592 fBgrIn = new TMap();
593 fBgrErrUncorrInSq = new TMap();
594 fBgrErrScaleIn = new TMap();
595 fSysIn = new TMap();
600 // results
601 fEmatUncorrX = 0;
602 fEmatUncorrAx = 0;
603 fDeltaCorrX = new TMap();
604 fDeltaCorrAx = new TMap();
607 fDeltaSysTau = 0;
608 fDtau=0.0;
609 fYData=0;
610 fVyyData=0;
611}
612
613////////////////////////////////////////////////////////////////////////////////
614/// Clear all data members which depend on the unfolding results.
615
617{
624}
625
626////////////////////////////////////////////////////////////////////////////////
627/// Matrix calculations required to propagate systematic errors.
628///
629/// data members modified:
630/// fEmatUncorrX, fEmatUncorrAx, fDeltaCorrX, fDeltaCorrAx
631
633{
634 if(!fEmatUncorrX) {
636 }
637 TMatrixDSparse *AM0=0,*AM1=0;
638 if(!fEmatUncorrAx) {
639 if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
640 if(!AM1) {
642 Int_t *rows_cols=new Int_t[GetNy()];
643 Double_t *data=new Double_t[GetNy()];
644 for(Int_t i=0;i<GetNy();i++) {
645 rows_cols[i]=i;
646 data[i]=1.0;
647 }
649 (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
650 delete[] data;
651 delete[] rows_cols;
652 AddMSparse(AM1,-1.,one);
653 DeleteMatrix(&one);
655 }
656 }
657 if((!fDeltaSysTau )&&(fDtau>0.0)) {
662 for(Int_t i=0;i<n;i++) {
663 data[i] *= scale;
664 }
665 }
666
667 TMapIter sysErrIn(fSysIn);
668 const TObjString *key;
669
670 // calculate individual systematic errors
671 for(key=(const TObjString *)sysErrIn.Next();key;
672 key=(const TObjString *)sysErrIn.Next()) {
673 const TMatrixDSparse *dsys=
674 (const TMatrixDSparse *)((const TPair *)*sysErrIn)->Value();
675 const TPair *named_emat=(const TPair *)
677 if(!named_emat) {
679 fDeltaCorrX->Add(new TObjString(*key),emat);
680 }
681 named_emat=(const TPair *)fDeltaCorrAx->FindObject(key->GetString());
682 if(!named_emat) {
683 if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
684 if(!AM1) {
686 Int_t *rows_cols=new Int_t[GetNy()];
687 Double_t *data=new Double_t[GetNy()];
688 for(Int_t i=0;i<GetNy();i++) {
689 rows_cols[i]=i;
690 data[i]=1.0;
691 }
693 (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
694 delete[] data;
695 delete[] rows_cols;
696 AddMSparse(AM1,-1.,one);
697 DeleteMatrix(&one);
699 }
700 TMatrixDSparse *emat=PrepareCorrEmat(AM0,AM1,dsys);
701 fDeltaCorrAx->Add(new TObjString(*key),emat);
702 }
703 }
704 DeleteMatrix(&AM0);
705 DeleteMatrix(&AM1);
706}
707
708////////////////////////////////////////////////////////////////////////////////
709/// Covariance contribution from uncorrelated uncertainties of the
710/// response matrix.
711///
712/// \param[inout] ematrix covariance matrix histogram
713/// \param[in] binMap mapping of histogram bins
714/// \param[in] clearEmat if true, ematrix is cleared prior to adding
715/// this covariance matrix contribution
716///
717/// This method propagates the uncertainties of the response matrix
718/// histogram, specified with the constructor, to the unfolding
719/// result. It is assumed that the entries of that histogram are
720/// bin-to-bin uncorrelated. In many cases this corresponds to the
721/// "Monte Carlo statistical uncertainties".
722///
723/// The array <b>binMap</b> is explained with the method GetOutput().
724/// The flag <b>clearEmat</b> may be used to add covariance matrices from
725/// several uncertainty sources.
726///
727/// data members modified:
728/// fVYAx, fESparse, fEAtV, fErrorAStat
729
731(TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
732{
734 ErrorMatrixToHist(ematrix,fEmatUncorrX,binMap,clearEmat);
735}
736
737////////////////////////////////////////////////////////////////////////////////
738/// Propagate uncorrelated systematic errors to a covariance matrix.
739///
740/// \param[in] m_0 coefficients for error propagation
741/// \param[in] m_1 coefficients for error propagation
742///
743/// Returns the covariance matrix, propagates uncorrelated systematic errors to
744/// a covariance matrix. m_0,m_1 are the coefficients (matrices) for propagating
745/// the errors.
746///
747/// The error matrix is calculated by standard error propagation, where the
748/// derivative of the result vector X wrt the matrix A is given by:
749///
750/// \f[ \frac{dX_k}{dA_{ij}} = M0_{kj} Z0_i - M1_{ki} Z1_j \f]
751///
752/// where:
753//
754/// the matrices M0 and M1 are arguments to this function
755/// the vectors Z0, Z1 : GetDXDAZ()
756///
757/// The matrix A is calculated from a matrix B as
758///
759/// \f[ A_{ij} = \frac{B_{ij}}{\sum_k B_{kj}} \f]
760///
761/// where k runs over additional indices of B, not present in A.
762/// (underflow and overflow bins, used for efficiency corrections)
763///
764/// define: \f$ Norm_j = \sum_k B_{kj} \f$ (data member fSumOverY)
765///
766/// the derivative of A wrt this input matrix B is given by:
767///
768/// \f[ \frac{dA_{ij}}{dB_{kj}} = (\delta_{ik} - A_{ij} ) \frac{1}{Norm_j} \f]
769///
770/// The covariance matrix Vxx is:
771///
772/// \f[ Vxx_{mn} = \sum_{ijlk} \big[ (\frac{dX_m}{dA_{ij}}) (\frac{dA_{ij}}{dB_{}kj}) DB_{kj} (\frac{dX_n}{dA_{lj}}) (\frac{dA_{lj}}{dB_{kj}}) \big] \f]
773///
774/// where \f$ DB_{kj} \f$ is the error on \f$ B_{kj} \f$ squared.
775///
776/// Simplify the sum over k:
777///
778/// \f[ \sum_k \big[ (\frac{dA_{ij}}{dB_{kj}}) DB_{kj} (\frac{dA_{lj}}{dB_{kj}}) \big]
779/// = \sum_k \big[ (\delta_{ik} - A_{ij} ) \frac{1}{Norm_j} DB_{kj} (\delta_{lk} - A_{lj} ) \frac{1}{Norm_j} \big]
780/// = \sum_k \big[ (\delta_{ik} \delta_{lk} - \delta_{ik} A_{lj} - \delta_{lk} A_{ij} + A_{ij} A_{lj} ) \frac{DB_{kj}}{Norm_j^2} \big] \f]
781///
782/// introduce normalized errors: \f$ Rsq_{kj} = \frac{DB_{kj}}{Norm_j^2} \f$
783///
784/// after summing over k:
785/// \f[ \delta_{ik} \delta_{lk} Rsq_{kj} \to \delta_{il} Rsq_{ij} \f]
786/// \f[ \delta_{ik} A_{lj} Rsq_{kj} \to A_{lj} Rsq_{ij} \f]
787/// \f[ \delta_{lk} A_{ij} Rsq_{kj} \to A_{ij} Rsq_{lj} \f]
788/// \f[ A_{ij} A_{lj} Rsq_{kj} \to A_{ij} A_{lj} \sum_k(Rsq_{kj}) \f]
789///
790/// introduce sum of normalized errors squared: \f$ SRsq_j = \sum_k(Rsq_{kj}) \f$
791///
792/// Note: \f$ Rsq_{ij} \f$ is stored as `fDAinRelSq` (excludes extra indices of B)
793/// and \f$ SRsq_j \f$ is stored as `fDAinColRelSq` (sum includes all indices of B)
794///
795/// \f[ Vxx_{nm} = \sum_{ijl} \big[ (\frac{dX_m}{dA_{ij}}) (\frac{dX_n}{dA_{lj}})
796/// (\delta_{il} Rsq_{ij} - A_{lj} Rsq_{ij} - A_{ij} Rsq_{lj} + A_{ij} A_{lj} SRsq_j) \big] \f]
797///
798/// \f[ Vxx_nm = \sum_j \big[ F_{mj} F_{nj} SRsq_j \big]
799/// - \sum_j \big[ G_{mj} F_{nj} \big]
800/// - \sum_j \big[ F_{mj} G_{nj} \big]
801/// + \sum_{ij} \big[ (\frac{dX_m}{dA_{ij}}) (\frac{dX_n}{dA_{lj}}) Rsq_{ij} \big] \f]
802///
803/// where:
804///
805/// \f[ F_{mj} = \sum_i \big[ (\frac{dX_m}{dA_{ij}}) * A_{ij} \big] \f]
806/// \f[ G_{mj} = \sum_i \big[ (\frac{dX_m}{dA_{ij}}) Rsq_{ij} \big] \f]
807///
808/// In order to avoid explicitly calculating the 3-dimensional tensor
809/// \f$(\frac{dX_m}{dA_{ij}}) \f$ the sums are evaluated further, using:
810///
811/// \f[ \frac{dX_k}{dA_{ij}} = M0_{kj} Z0_i - M1_{ki} Z1_j \f]
812/// \f[ F_{mj} = M0_{mj} * (A\# Z0)_j - (M1 A)_{mj} Z1_j \f]
813/// \f[ G_{mj} = M0_{mj} * (Rsq\# Z0)_j - (M1 Rsq)_{mj} Z1_j \f]
814///
815/// and
816///
817/// \f[ \sum_{ij} \big[ (\frac{dX_m}{dA_{ij}}) (\frac{dX_n}{dA_{ij}}) Rsq_{ij} \big] =
818/// \sum_j \big[ M0_{mj} M0_nj \big[ \sum_i (Z0_i)^2 Rsq_{ij} \big] \big]
819/// + \sum_i \big[ M1_{mi} M1_{ni} \big[ \sum_j (Z1_j)^2 Rsq_{ij} \big] \big]
820/// - \sum_i \big[ M1_{mi} H_{ni} + M1_{ni} H_{mi} \big] \f]
821///
822/// where:
823///
824/// \f[ H_{mi} = Z0_i \sum_j \big[ M0_{mj} Z1_j Rsq_{ij} \big] \f]
825///
826/// collect all contributions:
827///
828/// \f[ Vxx_nm = r0 -r1 -r2 +r3 +r4 -r5 -r6 \f]
829/// \f[ r0 = \sum_j \big[ F_{mj} F_nj * SRsq_j \big] \f]
830/// \f[ r1 = \sum_j \big[ G_{mj} F_nj \big] \f]
831/// \f[ r2 = \sum_j \big[ F_{mj} G_nj \big] \f]
832/// \f[ r3 = \sum_j \big[ M0_{mj} M0_nj \big[ \sum_i (Z0_i)^2 Rsq_{ij} \big] \big] \f]
833/// \f[ r4 = \sum_i \big[ M1_{mi} M1_{ni} \big[ \sum_j (Z1_j)^2 Rsq_{ij} \big] \big] \f]
834/// \f[ r5 = \sum_i \big[ M1_{mi} H_{ni} \big] \f]
835/// \f[ r6 = \sum_i \big[ M1_{ni} H_{mi} \big] \f]
836
838(const TMatrixDSparse *m_0,const TMatrixDSparse *m_1)
839{
840
841 //======================================================
842 // calculate contributions containing matrices F and G
843 // r0,r1,r2
844 TMatrixDSparse *r=0;
846 // calculate matrices (M1*A)_{mj} * Z1_j and (M1*Rsq)_{mj} * Z1_j
850 ScaleColumnsByVector(M1Rsq_Z1,GetDXDAZ(1));
851 // calculate vectors A#*Z0 and Rsq#*Z0
853 TMatrixDSparse *RsqZ0=
855 //calculate matrix F
856 // F_{mj} = M0_{mj} * (A# Z0)_j - (M1 A)_{mj} Z1_j
859 AddMSparse(F,-1.0,M1A_Z1);
860 //calculate matrix G
861 // G_{mj} = M0_{mj} * (Rsq# Z0)_j - (M1 Rsq)_{mj} Z1_j
863 ScaleColumnsByVector(G,RsqZ0);
864 AddMSparse(G,-1.0,M1Rsq_Z1);
865 DeleteMatrix(&M1A_Z1);
866 DeleteMatrix(&M1Rsq_Z1);
867 DeleteMatrix(&AtZ0);
868 DeleteMatrix(&RsqZ0);
869 // r0 = \sum_j [ F_{mj} * F_nj * SRsq_j ]
871 // r1 = \sum_j [ G_{mj} * F_nj ]
873 // r2 = \sum_j [ F_{mj} * G_nj ]
875 // r = r0-r1-r2
876 AddMSparse(r,-1.0,r1);
877 AddMSparse(r,-1.0,r2);
878 DeleteMatrix(&r1);
879 DeleteMatrix(&r2);
880 DeleteMatrix(&F);
881 DeleteMatrix(&G);
882 }
883 //======================================================
884 // calculate contribution
885 // \sum_{ij} [ (dX_m/dA_{ij}) * (dX_n/dA_{ij}) * Rsq_{ij} ]
886 // (r3,r4,r5,r6)
887 if(fDAinRelSq) {
888 // (Z0_i)^2
889 TMatrixDSparse Z0sq(*GetDXDAZ(0));
890 const Int_t *Z0sq_rows=Z0sq.GetRowIndexArray();
891 Double_t *Z0sq_data=Z0sq.GetMatrixArray();
892 for(int index=0;index<Z0sq_rows[Z0sq.GetNrows()];index++) {
893 Z0sq_data[index] *= Z0sq_data[index];
894 }
895 // Z0sqRsq = \sum_i (Z_i)^2 * Rsq_{ij}
897 // r3 = \sum_j [ M0_{mj} * M0_nj * [ \sum_i (Z0_i)^2 * Rsq_{ij} ] ]
899 DeleteMatrix(&Z0sqRsq);
900
901 // (Z1_j)^2
902 TMatrixDSparse Z1sq(*GetDXDAZ(1));
903 const Int_t *Z1sq_rows=Z1sq.GetRowIndexArray();
904 Double_t *Z1sq_data=Z1sq.GetMatrixArray();
905 for(int index=0;index<Z1sq_rows[Z1sq.GetNrows()];index++) {
906 Z1sq_data[index] *= Z1sq_data[index];
907 }
908 // Z1sqRsq = \sum_j (Z1_j)^2 * Rsq_{ij} ]
910 // r4 = \sum_i [ M1_{mi} * M1_{ni} * [ \sum_j (Z1_j)^2 * Rsq_{ij} ] ]
912 DeleteMatrix(&Z1sqRsq);
913
914 // \sum_j [ M0_{mj} * Z1_j * Rsq_{ij} ]
916 (m_0,fDAinRelSq,GetDXDAZ(1));
917 // H_{mi} = Z0_i * \sum_j [ M0_{mj} * Z1_j * Rsq_{ij} ]
919 // r5 = \sum_i [ M1_{mi} * H_{ni} ]
921 // r6 = \sum_i [ H_{mi} * M1_{ni} ]
923 DeleteMatrix(&H);
924 // r = r0 -r1 -r2 +r3 +r4 -r5 -r6
925 if(r) {
926 AddMSparse(r,1.0,r3);
927 DeleteMatrix(&r3);
928 } else {
929 r=r3;
930 r3=0;
931 }
932 AddMSparse(r,1.0,r4);
933 AddMSparse(r,-1.0,r5);
934 AddMSparse(r,-1.0,r6);
935 DeleteMatrix(&r4);
936 DeleteMatrix(&r5);
937 DeleteMatrix(&r6);
938 }
939 return r;
940}
941
942////////////////////////////////////////////////////////////////////////////////
943/// Propagate correlated systematic shift to an output vector.
944///
945/// \param[in] m1 coefficients
946/// \param[in] m2 coeffiicients
947/// \param[in] dsys matrix of correlated shifts from this source
948/// propagate correlated systematic shift to output vector
949/// m1,m2 : coefficients for propagating the errors
950/// dsys : matrix of correlated shifts from this source
951///
952/// \f[ \delta_m =
953/// \sum{i,j} {
954/// ((*m1)(m,j) * (*fVYAx)(i) - (*m2)(m,i) * (*fX)(j))*dsys(i,j) }
955/// = \sum_j (*m1)(m,j) \sum_i dsys(i,j) * (*fVYAx)(i)
956/// - \sum_i (*m2)(m,i) \sum_j dsys(i,j) * (*fX)(j) \f]
957
959(const TMatrixDSparse *m1,const TMatrixDSparse *m2,const TMatrixDSparse *dsys)
960{
962 TMatrixDSparse *delta = MultiplyMSparseMSparse(m1,dsysT_VYAx);
963 DeleteMatrix(&dsysT_VYAx);
965 TMatrixDSparse *delta2 = MultiplyMSparseMSparse(m2,dsys_X);
966 DeleteMatrix(&dsys_X);
967 AddMSparse(delta,-1.0,delta2);
968 DeleteMatrix(&delta2);
969 return delta;
970}
971
972////////////////////////////////////////////////////////////////////////////////
973/// Specify an uncertainty on tau.
974///
975/// \param[in] delta_tau new uncertainty on tau
976///
977/// The default is to have no uncertyainty on tau.
978
980{
981 fDtau=delta_tau;
983}
984
985////////////////////////////////////////////////////////////////////////////////
986/// Correlated one-sigma shifts correspinding to a given systematic uncertainty.
987///
988/// \param[out] hist_delta histogram to store shifts
989/// \param[in] name identifier of the background source
990/// \param[in] binMap (default=0) remapping of histogram bins
991///
992/// returns true if the error source was found.
993///
994/// This method returns the shifts of the unfolding result induced by
995/// varying the identified systematic source by one sigma.
996///
997/// the array <b>binMap</b> is explained with the method GetOutput().
998
1000 const Int_t *binMap)
1001{
1003 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1004 const TMatrixDSparse *delta=0;
1005 if(named_emat) {
1006 delta=(TMatrixDSparse *)named_emat->Value();
1007 }
1008 VectorMapToHist(hist_delta,delta,binMap);
1009 return delta !=0;
1010}
1011
1012////////////////////////////////////////////////////////////////////////////////
1013/// Correlated one-sigma shifts from background normalisation uncertainty.
1014///
1015/// \param[out] hist_delta histogram to store shifts
1016/// \param[in] source identifier of the background source
1017/// \param[in] binMap (default=0) remapping of histogram bins
1018///
1019/// returns true if the background source was found.
1020///
1021/// This method returns the shifts of the unfolding result induced by
1022/// varying the normalisation of the identified background by one sigma.
1023///
1024/// the array <b>binMap</b> is explained with the method GetOutput().
1025
1027(TH1 *hist_delta,const char *source,const Int_t *binMap)
1028{
1030 const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(source);
1031 TMatrixDSparse *dx=0;
1032 if(named_err) {
1033 const TMatrixD *dy=(TMatrixD *)named_err->Value();
1034 dx=MultiplyMSparseM(GetDXDY(),dy);
1035 }
1036 VectorMapToHist(hist_delta,dx,binMap);
1037 if(dx!=0) {
1038 DeleteMatrix(&dx);
1039 return kTRUE;
1040 }
1041 return kFALSE;
1042}
1043
1044////////////////////////////////////////////////////////////////////////////////
1045/// Correlated one-sigma shifts from shifting tau.
1046///
1047/// \param[out] hist_delta histogram to store shifts
1048/// \param[in] binMap (default=0) remapping of histogram bins
1049///
1050/// returns true if the background source was found.
1051///
1052/// This method returns the shifts of the unfolding result induced by
1053/// varying the normalisation of the identified background by one sigma.
1054///
1055/// the array <b>binMap</b> is explained with the method GetOutput().
1056///
1057/// calculate systematic shift from tau variation
1058/// - ematrix: output
1059/// - binMap: see method GetEmatrix()
1060
1061Bool_t TUnfoldSys::GetDeltaSysTau(TH1 *hist_delta,const Int_t *binMap)
1062{
1064 VectorMapToHist(hist_delta,fDeltaSysTau,binMap);
1065 return fDeltaSysTau !=0;
1066}
1067
1068////////////////////////////////////////////////////////////////////////////////
1069/// Covariance contribution from a systematic variation of the
1070/// response matrix.
1071///
1072/// \param[inout] ematrix covariance matrix histogram
1073/// \param[in] name identifier of the systematic variation
1074/// \param[in] binMap (default=0) remapping of histogram bins
1075/// \param[in] clearEmat (default=true) if true, clear the histogram
1076/// prior to adding the covariance matrix contribution
1077///
1078/// Returns the covariance matrix contribution from shifting the given
1079/// uncertainty source within one sigma
1080///
1081/// the array <b>binMap</b> is explained with the method GetOutput().
1082/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1083/// several uncertainty sources.
1084
1086(TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1087{
1089 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1090 TMatrixDSparse *emat=0;
1091 if(named_emat) {
1092 const TMatrixDSparse *delta=(TMatrixDSparse *)named_emat->Value();
1093 emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
1094 }
1095 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1096 DeleteMatrix(&emat);
1097}
1098
1099////////////////////////////////////////////////////////////////////////////////
1100/// Covariance contribution from background normalisation uncertainty.
1101///
1102/// \param[in,out] ematrix output histogram
1103/// \param[in] name identifier of the background source
1104/// \param[in] binMap (default=0) remapping of histogram bins
1105/// \param[in] clearEmat (default=true) if true, clear the histogram
1106/// prior to adding the covariance matrix contribution
1107///
1108/// this method returns the uncertainties on the unfolding result
1109/// arising from the background source <b>source</b> and its normalisation
1110/// uncertainty. See method SubtractBackground() how to set the normalisation uncertainty
1111///
1112/// the array <b>binMap</b> is explained with the method GetOutput().
1113/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1114/// several uncertainty sources.
1115
1117(TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1118{
1120 const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(name);
1121 TMatrixDSparse *emat=0;
1122 if(named_err) {
1123 const TMatrixD *dy=(TMatrixD *)named_err->Value();
1126 DeleteMatrix(&dx);
1127 }
1128 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1129 DeleteMatrix(&emat);
1130}
1131
1132////////////////////////////////////////////////////////////////////////////////
1133/// Covariance matrix contribution from error on regularisation
1134/// parameter.
1135///
1136/// \param[inout] ematrix output histogram
1137/// \param[in] binMap (default=0) remapping of histogram bins
1138/// \param[in] clearEmat (default=true) if true, clear the histogram
1139///
1140/// this method returns the covariance contributions to the unfolding result
1141/// from the assigned uncertainty on the parameter tau, see method
1142/// SetTauError().
1143///
1144/// the array <b>binMap</b> is explained with the method GetOutput().
1145/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1146/// several uncertainty sources.
1147///
1148/// Calculate error matrix from error in regularisation parameter
1149/// - ematrix: output
1150/// - binMap: see method GetEmatrix()
1151/// - clearEmat: set kTRUE to clear the histogram prior to adding the errors
1152
1154(TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1155{
1157 TMatrixDSparse *emat=0;
1158 if(fDeltaSysTau) {
1160 }
1161 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1162 DeleteMatrix(&emat);
1163}
1164
1165////////////////////////////////////////////////////////////////////////////////
1166/// Covariance matrix contribution from input measurement uncertainties.
1167///
1168/// \param[inout] ematrix output histogram
1169/// \param[in] binMap (default=0) remapping of histogram bins
1170/// \param[in] clearEmat (default=true) if true, clear the histogram
1171///
1172/// this method returns the covariance contributions to the unfolding result
1173/// from the uncertainties or covariance of the input
1174/// data. In many cases, these are the "statistical uncertainties".
1175///
1176/// The array <b>binMap</b> is explained with the method GetOutput().
1177/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1178/// several uncertainty sources.
1179
1181(TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1182{
1183 GetEmatrixFromVyy(fVyyData,ematrix,binMap,clearEmat);
1184}
1185
1186////////////////////////////////////////////////////////////////////////////////
1187/// Covariance contribution from background uncorrelated uncertainty.
1188///
1189/// \param[in] ematrix output histogram
1190/// \param[in] source identifier of the background source
1191/// \param[in] binMap (default=0) remapping of histogram bins
1192/// \param[in] clearEmat (default=true) if true, clear the histogram
1193///
1194/// this method returns the covariance contributions to the unfolding result
1195/// arising from the background source <b>source</b> and the uncorrelated
1196/// (background histogram uncertainties). Also see method SubtractBackground()
1197///
1198/// the array <b>binMap</b> is explained with the method GetOutput().
1199/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1200/// several uncertainty sources.
1201
1203(TH2 *ematrix,const char *source,const Int_t *binMap,Bool_t clearEmat)
1204{
1205 const TPair *named_err=(const TPair *)fBgrErrUncorrInSq->FindObject(source);
1206 TMatrixDSparse *emat=0;
1207 if(named_err) {
1208 TMatrixD const *dySquared=(TMatrixD const *)named_err->Value();
1210 }
1211 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1212 DeleteMatrix(&emat);
1213}
1214
1215////////////////////////////////////////////////////////////////////////////////
1216/// Propagate an error matrix on the input vector to the unfolding result.
1217///
1218/// \param[in] vyy input error matrix
1219/// \param[inout] ematrix histogram to be updated
1220/// \param[in] binMap mapping of histogram bins
1221/// \param[in] clearEmat if set, clear histogram before adding this
1222/// covariance contribution
1223///
1224/// propagate error matrix vyy to the result
1225/// - vyy: error matrix on input data fY
1226/// - ematrix: output
1227/// - binMap: see method GetEmatrix()
1228/// - clearEmat: set kTRUE to clear the histogram prior to adding the errors
1229
1231(const TMatrixDSparse *vyy,TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1232{
1234 TMatrixDSparse *em=0;
1235 if(vyy) {
1238 DeleteMatrix(&dxdyVyy);
1239 }
1240 ErrorMatrixToHist(ematrix,em,binMap,clearEmat);
1241 DeleteMatrix(&em);
1242}
1243
1244////////////////////////////////////////////////////////////////////////////////
1245/// Get total error matrix, summing up all contributions.
1246///
1247/// \param[out] ematrix histogram which will be filled
1248/// \param[in] binMap (default=0) remapping of histogram bins
1249///
1250/// the array <b>binMap</b> is explained with the method GetOutput().
1251///
1252/// get total error including statistical error
1253/// - ematrix: output
1254/// - binMap: see method GetEmatrix()
1255
1256void TUnfoldSys::GetEmatrixTotal(TH2 *ematrix,const Int_t *binMap)
1257{
1258 GetEmatrix(ematrix,binMap); // (stat)+(d)+(e)
1259 GetEmatrixSysUncorr(ematrix,binMap,kFALSE); // (a)
1260 TMapIter sysErrPtr(fDeltaCorrX);
1261 const TObject *key;
1262
1263 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1264 GetEmatrixSysSource(ematrix,
1265 ((const TObjString *)key)->GetString(),
1266 binMap,kFALSE); // (b)
1267 }
1268 GetEmatrixSysTau(ematrix,binMap,kFALSE); // (c)
1269}
1270
1271////////////////////////////////////////////////////////////////////////////////
1272/// Determine total error matrix on the vector Ax.
1273
1275{
1277
1278 // errors from input vector and from background subtraction
1279 TMatrixDSparse *emat_sum=new TMatrixDSparse(*fVyy);
1280
1281 // uncorrelated systematic error
1282 if(fEmatUncorrAx) {
1283 AddMSparse(emat_sum,1.0,fEmatUncorrAx);
1284 }
1285 TMapIter sysErrPtr(fDeltaCorrAx);
1286 const TObject *key;
1287
1288 // correlated systematic errors
1289 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1290 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1292 AddMSparse(emat_sum,1.0,emat);
1293 DeleteMatrix(&emat);
1294 }
1295 // error on tau
1296 if(fDeltaSysTau) {
1298 TMatrixDSparse *emat_tau=
1299 MultiplyMSparseMSparseTranspVector(Adx_tau,Adx_tau,0);
1300 DeleteMatrix(&Adx_tau);
1301 AddMSparse(emat_sum,1.0,emat_tau);
1302 DeleteMatrix(&emat_tau);
1303 }
1304 return emat_sum;
1305}
1306
1307////////////////////////////////////////////////////////////////////////////////
1308/// Determine total error matrix on the vector x.
1309
1311{
1313
1314 // errors from input vector and from background subtraction
1315 TMatrixDSparse *emat_sum=new TMatrixDSparse(*GetVxx());
1316
1317 // uncorrelated systematic error
1318 if(fEmatUncorrX) {
1319 AddMSparse(emat_sum,1.0,fEmatUncorrX);
1320 }
1321 TMapIter sysErrPtr(fDeltaCorrX);
1322 const TObject *key;
1323
1324 // correlated systematic errors
1325 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1326 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1328 AddMSparse(emat_sum,1.0,emat);
1329 DeleteMatrix(&emat);
1330 }
1331 // error on tau
1332 if(fDeltaSysTau) {
1333 TMatrixDSparse *emat_tau=
1335 AddMSparse(emat_sum,1.0,emat_tau);
1336 DeleteMatrix(&emat_tau);
1337 }
1338 return emat_sum;
1339}
1340
1341
1342////////////////////////////////////////////////////////////////////////////////
1343/// Calculate total chi**2 including all systematic errors.
1344
1346{
1347
1349
1350 Int_t rank=0;
1351 TMatrixDSparse *v=InvertMSparseSymmPos(emat_sum,&rank);
1354 DeleteMatrix(&v);
1355 Double_t r=0.0;
1356 const Int_t *vdy_rows=vdy->GetRowIndexArray();
1357 const Double_t *vdy_data=vdy->GetMatrixArray();
1358 for(Int_t i=0;i<vdy->GetNrows();i++) {
1359 if(vdy_rows[i+1]>vdy_rows[i]) {
1360 r += vdy_data[vdy_rows[i]] * dy(i,0);
1361 }
1362 }
1363 DeleteMatrix(&vdy);
1364 DeleteMatrix(&emat_sum);
1365 return r;
1366}
1367
1368////////////////////////////////////////////////////////////////////////////////
1369/// Get global correlatiocn coefficients, summing up all contributions.
1370///
1371/// \param[out] rhoi histogram which will be filled
1372/// \param[in] binMap (default=0) remapping of histogram bins
1373/// \param[out] invEmat (default=0) inverse of error matrix
1374///
1375/// return the global correlation coefficients, including all error
1376/// sources. If <b>invEmat</b> is nonzero, the inverse of the error
1377/// matrix is returned in that histogram
1378///
1379/// the array <b>binMap</b> is explained with the method GetOutput().
1380///
1381/// get global correlation coefficients including systematic,statistical,background,tau errors
1382/// - rhoi: output histogram
1383/// - binMap: for each global bin, indicate in which histogram bin
1384/// to store its content
1385/// - invEmat: output histogram for inverse of error matrix
1386/// (pointer may zero if inverse is not requested)
1387
1388void TUnfoldSys::GetRhoItotal(TH1 *rhoi,const Int_t *binMap,TH2 *invEmat)
1389{
1390 ClearHistogram(rhoi,-1.);
1392 GetRhoIFromMatrix(rhoi,emat_sum,binMap,invEmat);
1393
1394 DeleteMatrix(&emat_sum);
1395}
1396
1397////////////////////////////////////////////////////////////////////////////////
1398/// Scale columns of a matrix by the corresponding rows of a vector.
1399///
1400/// \param[inout] m matrix
1401/// \param[in] v vector
1402///
1403/// the entries m<sub>ij</sub> are multiplied by v<sub>j</sub>.
1404///
1405/// scale columns of m by the corresponding rows of v
1406/// input:
1407/// - m: pointer to sparse matrix of dimension NxM
1408/// - v: pointer to matrix of dimension Mx1
1409
1412{
1413 if((m->GetNcols() != v->GetNrows())||(v->GetNcols()!=1)) {
1414 Fatal("ScaleColumnsByVector error",
1415 "matrix cols/vector rows %d!=%d OR vector cols %d !=1\n",
1416 m->GetNcols(),v->GetNrows(),v->GetNcols());
1417 }
1418 const Int_t *rows_m=m->GetRowIndexArray();
1419 const Int_t *cols_m=m->GetColIndexArray();
1420 Double_t *data_m=m->GetMatrixArray();
1421 const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
1422 if(v_sparse) {
1423 const Int_t *rows_v=v_sparse->GetRowIndexArray();
1424 const Double_t *data_v=v_sparse->GetMatrixArray();
1425 for(Int_t i=0;i<m->GetNrows();i++) {
1426 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1427 Int_t j=cols_m[index_m];
1428 Int_t index_v=rows_v[j];
1429 if(index_v<rows_v[j+1]) {
1430 data_m[index_m] *= data_v[index_v];
1431 } else {
1432 data_m[index_m] =0.0;
1433 }
1434 }
1435 }
1436 } else {
1437 for(Int_t i=0;i<m->GetNrows();i++) {
1438 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1439 Int_t j=cols_m[index_m];
1440 data_m[index_m] *= (*v)(j,0);
1441 }
1442 }
1443 }
1444}
1445
1446////////////////////////////////////////////////////////////////////////////////
1447/// Map delta to hist_delta, possibly summing up bins.
1448///
1449/// \param[out] hist_delta result histogram
1450/// \param[in] delta vector to be mapped to the histogram
1451/// \param[in] binMap mapping of histogram bins
1452///
1453/// groups of bins of <b>delta</b> are mapped to bins of
1454/// <b>hist_delta</b>. The histogram contents are set to the sum over
1455/// the group of bins. The histogram errors are reset to zero.
1456///
1457/// The array <b>binMap</b> is explained with the method GetOutput()
1458///
1459/// sum over bins of *delta, as defined in binMap,fXToHist
1460/// - hist_delta: histogram to return summed vector
1461/// - delta: vector to sum and remap
1462
1464(TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap)
1465{
1466 Int_t nbin=hist_delta->GetNbinsX();
1467 Double_t *c=new Double_t[nbin+2];
1468 for(Int_t i=0;i<nbin+2;i++) {
1469 c[i]=0.0;
1470 }
1471 if(delta) {
1472 Int_t binMapSize = fHistToX.GetSize();
1473 const Double_t *delta_data=delta->GetMatrixArray();
1474 const Int_t *delta_rows=delta->GetRowIndexArray();
1475 for(Int_t i=0;i<binMapSize;i++) {
1476 Int_t destBinI=binMap ? binMap[i] : i;
1477 Int_t srcBinI=fHistToX[i];
1478 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
1479 Int_t index=delta_rows[srcBinI];
1480 if(index<delta_rows[srcBinI+1]) {
1481 c[destBinI]+=delta_data[index];
1482 }
1483 }
1484 }
1485 }
1486 for(Int_t i=0;i<nbin+2;i++) {
1487 hist_delta->SetBinContent(i,c[i]);
1488 hist_delta->SetBinError(i,0.0);
1489 }
1490 delete[] c;
1491}
1492
1493////////////////////////////////////////////////////////////////////////////////
1494/// Get a new list of all systematic uuncertainty sources.
1495///
1496/// The user is responsible for deleting the list
1497/// get list of names of systematic sources
1498
1500 TSortedList *r=new TSortedList();
1501 TMapIter i(fSysIn);
1502 for(const TObject *key=i.Next();key;key=i.Next()) {
1503 r->Add(((TObjString *)key)->Clone());
1504 }
1505 return r;
1506}
1507
1508////////////////////////////////////////////////////////////////////////////////
1509/// Get a new list of all background sources.
1510///
1511/// The user is responsible for deleting the list
1512/// get list of name of background sources
1513
1515 TSortedList *r=new TSortedList();
1516 TMapIter i(fBgrIn);
1517 for(const TObject *key=i.Next();key;key=i.Next()) {
1518 r->Add(((TObjString *)key)->Clone());
1519 }
1520 return r;
1521}
#define c(i)
Definition: RSha256.hxx:101
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
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
Definition: TCollection.h:179
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8929
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
virtual void SetBinError(Int_t bin, Double_t error)
Set the bin Error Note that this resets the bin eror option to be of Normal Type and for the non-empt...
Definition: TH1.cxx:9072
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:9088
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
Service class for 2-D histogram classes.
Definition: TH2.h:30
Double_t GetBinContent(Int_t binx, Int_t biny) const override
Definition: TH2.h:89
Iterator of map.
Definition: TMap.h:144
TObject * Next() override
Returns the next key from a map.
Definition: TMap.cxx:542
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:54
virtual void SetOwnerKeyValue(Bool_t ownkeys=kTRUE, Bool_t ownvals=kTRUE)
Set ownership for keys and values.
Definition: TMap.cxx:352
TObject * FindObject(const char *keyname) const override
Check if a (key,value) pair exists with keyname as name of the key.
Definition: TMap.cxx:215
void Clear(Option_t *option="") override
Remove all (key,value) pairs from the map.
Definition: TMap.cxx:97
Int_t GetNrows() const
Definition: TMatrixTBase.h:123
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 TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:223
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:955
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:969
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:997
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:943
Class used by TMap to store (key,value) pairs.
Definition: TMap.h:102
TObject * Value() const
Definition: TMap.h:121
A sorted doubly linked list.
Definition: TSortedList.h:28
Basic string class.
Definition: TString.h:136
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:451
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
TMatrixD * fAoutside
Input: underflow/overflow bins.
Definition: TUnfoldSys.h:64
TMatrixDSparse * fDAinRelSq
Input: normalized errors from input matrix.
Definition: TUnfoldSys.h:60
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:86
TMatrixD * fYData
Input: fY prior to bgr subtraction.
Definition: TUnfoldSys.h:76
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.
Definition: TUnfoldSys.cxx:530
void InitTUnfoldSys(void)
Initialize pointers and TMaps.
Definition: TUnfoldSys.cxx:586
TMatrixDSparse * fVyyData
Input: error on fY prior to bgr subtraction.
Definition: TUnfoldSys.h:78
TMatrixDSparse * fEmatUncorrAx
Result: syst.error from fDA2 on fAx.
Definition: TUnfoldSys.h:82
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.
Definition: TUnfoldSys.cxx:333
TMatrixDSparse * GetSummedErrorMatrixYY(void)
Determine total error matrix on the vector Ax.
Double_t fDtau
Input: error on tau.
Definition: TUnfoldSys.h:74
ESysErrMode
type of matrix specified with AddSysError()
Definition: TUnfoldSys.h:102
@ kSysErrModeRelative
matrix gives the relative shifts
Definition: TUnfoldSys.h:108
@ kSysErrModeMatrix
matrix is an alternative to the default matrix, the errors are the difference to the original matrix
Definition: TUnfoldSys.h:104
@ kSysErrModeShift
matrix gives the absolute shifts
Definition: TUnfoldSys.h:106
virtual TMatrixDSparse * PrepareUncorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2)
Propagate uncorrelated systematic errors to a covariance matrix.
Definition: TUnfoldSys.cxx:838
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.
Definition: TUnfoldSys.cxx:999
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:80
~TUnfoldSys(void) override
Definition: TUnfoldSys.cxx:119
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.
Definition: TUnfoldSys.cxx:483
void ClearResults(void) override
Clear all data members which depend on the unfolding results.
Definition: TUnfoldSys.cxx:616
TMap * fBgrIn
Input: size of background sources.
Definition: TUnfoldSys.h:68
TMap * fDeltaCorrX
Result: syst.shift from fSysIn on fX.
Definition: TUnfoldSys.h:84
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:62
TMatrixDSparse * fDeltaSysTau
Result: systematic shift from tau.
Definition: TUnfoldSys.h:88
void AddSysError(const TH2 *sysError, const char *name, EHistMap histmap, ESysErrMode mode)
Specify a correlated systematic uncertainty.
Definition: TUnfoldSys.cxx:248
void SetTauError(Double_t delta_tau)
Specify an uncertainty on tau.
Definition: TUnfoldSys.cxx:979
TMap * fBgrErrUncorrInSq
Input: uncorr error squared from bgr sources.
Definition: TUnfoldSys.h:70
TUnfoldSys(void)
Only for use by root streamer or derived classes.
Definition: TUnfoldSys.cxx:138
virtual void PrepareSysError(void)
Matrix calculations required to propagate systematic errors.
Definition: TUnfoldSys.cxx:632
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 the input data for subsequent calls to DoUnfold(Double_t).
Definition: TUnfoldSys.cxx:444
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.
Definition: TUnfoldSys.cxx:731
TMap * fBgrErrScaleIn
Input: background sources correlated error.
Definition: TUnfoldSys.h:72
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.
Definition: TUnfoldSys.cxx:959
TMap * fSysIn
Input: correlated errors.
Definition: TUnfoldSys.h:66
An algorithm to unfold distributions from detector to truth level.
Definition: TUnfold.h:103
TArrayI fHistToX
mapping of histogram bins to matrix indices
Definition: TUnfold.h:166
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse *a, const TMatrixD *b) const
Multiply sparse matrix and a non-sparse matrix.
Definition: TUnfold.cxx:774
TMatrixDSparse * MultiplyMSparseTranspMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply a transposed Sparse matrix with another sparse matrix,.
Definition: TUnfold.cxx:694
TMatrixDSparse * MultiplyMSparseMSparseTranspVector(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixTBase< Double_t > *v) const
Calculate a sparse matrix product where the diagonal matrix V is given by a vector.
Definition: TUnfold.cxx:834
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:593
const TMatrixDSparse * GetDXDAM(int i) const
matrix contributions of the derivative dx/dA
Definition: TUnfold.h:248
Int_t GetNy(void) const
returns the number of measurement bins
Definition: TUnfold.h:234
const TMatrixDSparse * GetDXDtauSquared(void) const
vector of derivative dx/dtauSquared, using internal bin counting
Definition: TUnfold.h:261
static void DeleteMatrix(TMatrixD **m)
delete matrix and invalidate pointer
Definition: TUnfold.cxx:196
void ClearHistogram(TH1 *h, Double_t x=0.) const
Initialize bin contents and bin errors for a given histogram.
Definition: TUnfold.cxx:3644
Int_t GetNx(void) const
returns internal number of output (truth) matrix rows
Definition: TUnfold.h:226
const TMatrixDSparse * GetDXDAZ(int i) const
vector contributions of the derivative dx/dA
Definition: TUnfold.h:250
EConstraint
type of extra constraint
Definition: TUnfold.h:109
TMatrixDSparse * fVyy
covariance matrix Vyy corresponding to y
Definition: TUnfold.h:156
const TMatrixDSparse * GetVyyInv(void) const
inverse of covariance matrix of the data y
Definition: TUnfold.h:256
TArrayD fSumOverY
truth vector calculated from the non-normalized response matrix
Definition: TUnfold.h:168
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
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:3344
TArrayI fXToHist
mapping of matrix indices to histogram bins
Definition: TUnfold.h:164
const TMatrixDSparse * GetAx(void) const
vector of folded-back result
Definition: TUnfold.h:244
TMatrixD * fY
input (measured) data y
Definition: TUnfold.h:154
TMatrixDSparse * InvertMSparseSymmPos(const TMatrixDSparse *A, Int_t *rank) const
Get the inverse or pseudo-inverse of a positive, sparse matrix.
Definition: TUnfold.cxx:1009
Double_t fTauSquared
regularisation parameter tau squared
Definition: TUnfold.h:162
virtual void ClearResults(void)
Reset all results.
Definition: TUnfold.cxx:219
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=nullptr) const
Get output covariance matrix, possibly cumulated over several bins.
Definition: TUnfold.cxx:3411
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply two sparse matrices.
Definition: TUnfold.cxx:618
EHistMap
arrangement of axes for the response matrix (TH2 histogram)
Definition: TUnfold.h:139
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
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:931
const TMatrixDSparse * GetVxx(void) const
covariance matrix of the result
Definition: TUnfold.h:240
Double_t GetRhoIFromMatrix(TH1 *rhoi, const TMatrixDSparse *eOrig, const Int_t *binMap, TH2 *invEmat) const
Get global correlation coefficients with arbitrary min map.
Definition: TUnfold.cxx:3529
const TMatrixDSparse * GetDXDY(void) const
matrix of derivatives dx/dy
Definition: TUnfold.h:246
TMatrixDSparse * fA
response matrix A
Definition: TUnfold.h:150
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:2301
RVec< PromoteTypes< T0, T1 > > hypot(const T0 &x, const RVec< T1 > &v)
Definition: RVec.hxx:1795
const Int_t n
Definition: legend1.C:16
#define F(x, y, z)
#define G(x, y, z)
#define H(x, y, z)
static constexpr double m2
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:660
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
const char * Value
Definition: TXMLSetup.cxx:73
TMarker m
Definition: textangle.C:8
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345