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
278 if(mode!=kSysErrModeMatrix) {
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 }
287 if(mode==kSysErrModeShift) {
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] scaleError 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] source identifier of the background source
1049/// \param[in] binMap (default=0) remapping of histogram bins
1050///
1051/// returns true if the background source was found.
1052///
1053/// This method returns the shifts of the unfolding result induced by
1054/// varying the normalisation of the identified background by one sigma.
1055///
1056/// the array <b>binMap</b> is explained with the method GetOutput().
1057///
1058/// calculate systematic shift from tau variation
1059/// - ematrix: output
1060/// - binMap: see method GetEmatrix()
1061
1062Bool_t TUnfoldSys::GetDeltaSysTau(TH1 *hist_delta,const Int_t *binMap)
1063{
1065 VectorMapToHist(hist_delta,fDeltaSysTau,binMap);
1066 return fDeltaSysTau !=0;
1067}
1068
1069////////////////////////////////////////////////////////////////////////////////
1070/// Covariance contribution from a systematic variation of the
1071/// response matrix.
1072///
1073/// \param[inout] ematrix covariance matrix histogram
1074/// \param[in] name identifier of the systematic variation
1075/// \param[in] binMap (default=0) remapping of histogram bins
1076/// \param[in] clearEmat (default=true) if true, clear the histogram
1077/// prior to adding the covariance matrix contribution
1078///
1079/// Returns the covariance matrix contribution from shifting the given
1080/// uncertainty source within one sigma
1081///
1082/// the array <b>binMap</b> is explained with the method GetOutput().
1083/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1084/// several uncertainty sources.
1085
1087(TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1088{
1090 const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1091 TMatrixDSparse *emat=0;
1092 if(named_emat) {
1093 const TMatrixDSparse *delta=(TMatrixDSparse *)named_emat->Value();
1094 emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
1095 }
1096 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1097 DeleteMatrix(&emat);
1098}
1099
1100////////////////////////////////////////////////////////////////////////////////
1101/// Covariance contribution from background normalisation uncertainty.
1102///
1103/// \param[inout] ematrix output histogram
1104/// \param[in] source identifier of the background source
1105/// \param[in] binMap (default=0) remapping of histogram bins
1106/// \param[in] clearEmat (default=true) if true, clear the histogram
1107/// prior to adding the covariance matrix contribution
1108///
1109/// this method returns the uncertainties on the unfolding result
1110/// arising from the background source <b>source</b> and its normalisation
1111/// uncertainty. See method SubtractBackground() how to set the normalisation uncertainty
1112///
1113/// the array <b>binMap</b> is explained with the method GetOutput().
1114/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1115/// several uncertainty sources.
1116
1118(TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1119{
1121 const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(name);
1122 TMatrixDSparse *emat=0;
1123 if(named_err) {
1124 const TMatrixD *dy=(TMatrixD *)named_err->Value();
1127 DeleteMatrix(&dx);
1128 }
1129 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1130 DeleteMatrix(&emat);
1131}
1132
1133////////////////////////////////////////////////////////////////////////////////
1134/// Covariance matrix contribution from error on regularisation
1135/// parameter.
1136///
1137/// \param[inout] ematrix output histogram
1138/// \param[in] binMap (default=0) remapping of histogram bins
1139/// \param[in] clearEmat (default=true) if true, clear the histogram
1140///
1141/// this method returns the covariance contributions to the unfolding result
1142/// from the assigned uncertainty on the parameter tau, see method
1143/// SetTauError().
1144///
1145/// the array <b>binMap</b> is explained with the method GetOutput().
1146/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1147/// several uncertainty sources.
1148///
1149/// Calculate error matrix from error in regularisation parameter
1150/// - ematrix: output
1151/// - binMap: see method GetEmatrix()
1152/// - clearEmat: set kTRUE to clear the histogram prior to adding the errors
1153
1155(TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1156{
1158 TMatrixDSparse *emat=0;
1159 if(fDeltaSysTau) {
1161 }
1162 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1163 DeleteMatrix(&emat);
1164}
1165
1166////////////////////////////////////////////////////////////////////////////////
1167/// Covariance matrix contribution from input measurement uncertainties.
1168///
1169/// \param[inout] ematrix output histogram
1170/// \param[in] binMap (default=0) remapping of histogram bins
1171/// \param[in] clearEmat (default=true) if true, clear the histogram
1172///
1173/// this method returns the covariance contributions to the unfolding result
1174/// from the uncertainties or covariance of the input
1175/// data. In many cases, these are the "statistical uncertainties".
1176///
1177/// The array <b>binMap</b> is explained with the method GetOutput().
1178/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1179/// several uncertainty sources.
1180
1182(TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1183{
1184 GetEmatrixFromVyy(fVyyData,ematrix,binMap,clearEmat);
1185}
1186
1187////////////////////////////////////////////////////////////////////////////////
1188/// Covariance contribution from background uncorrelated uncertainty.
1189///
1190/// \param[in] ematrix output histogram
1191/// \param[in] source identifier of the background source
1192/// \param[in] binMap (default=0) remapping of histogram bins
1193/// \param[in] clearEmat (default=true) if true, clear the histogram
1194///
1195/// this method returns the covariance contributions to the unfolding result
1196/// arising from the background source <b>source</b> and the uncorrelated
1197/// (background histogram uncertainties). Also see method SubtractBackground()
1198///
1199/// the array <b>binMap</b> is explained with the method GetOutput().
1200/// The flag <b>clearEmat</b> may be used to add covariance matrices from
1201/// several uncertainty sources.
1202
1204(TH2 *ematrix,const char *source,const Int_t *binMap,Bool_t clearEmat)
1205{
1206 const TPair *named_err=(const TPair *)fBgrErrUncorrInSq->FindObject(source);
1207 TMatrixDSparse *emat=0;
1208 if(named_err) {
1209 TMatrixD const *dySquared=(TMatrixD const *)named_err->Value();
1211 }
1212 ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1213 DeleteMatrix(&emat);
1214}
1215
1216////////////////////////////////////////////////////////////////////////////////
1217/// Propagate an error matrix on the input vector to the unfolding result.
1218///
1219/// \param[in] vyy input error matrix
1220/// \param[inout] ematrix histogram to be updated
1221/// \param[in] binMap mapping of histogram bins
1222/// \param[in] clearEmat if set, clear histogram before adding this
1223/// covariance contribution
1224///
1225/// propagate error matrix vyy to the result
1226/// - vyy: error matrix on input data fY
1227/// - ematrix: output
1228/// - binMap: see method GetEmatrix()
1229/// - clearEmat: set kTRUE to clear the histogram prior to adding the errors
1230
1232(const TMatrixDSparse *vyy,TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1233{
1235 TMatrixDSparse *em=0;
1236 if(vyy) {
1239 DeleteMatrix(&dxdyVyy);
1240 }
1241 ErrorMatrixToHist(ematrix,em,binMap,clearEmat);
1242 DeleteMatrix(&em);
1243}
1244
1245////////////////////////////////////////////////////////////////////////////////
1246/// Get total error matrix, summing up all contributions.
1247///
1248/// \param[out] ematrix histogram which will be filled
1249/// \param[in] binMap (default=0) remapping of histogram bins
1250///
1251/// the array <b>binMap</b> is explained with the method GetOutput().
1252///
1253/// get total error including statistical error
1254/// - ematrix: output
1255/// - binMap: see method GetEmatrix()
1256
1257void TUnfoldSys::GetEmatrixTotal(TH2 *ematrix,const Int_t *binMap)
1258{
1259 GetEmatrix(ematrix,binMap); // (stat)+(d)+(e)
1260 GetEmatrixSysUncorr(ematrix,binMap,kFALSE); // (a)
1261 TMapIter sysErrPtr(fDeltaCorrX);
1262 const TObject *key;
1263
1264 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1265 GetEmatrixSysSource(ematrix,
1266 ((const TObjString *)key)->GetString(),
1267 binMap,kFALSE); // (b)
1268 }
1269 GetEmatrixSysTau(ematrix,binMap,kFALSE); // (c)
1270}
1271
1272////////////////////////////////////////////////////////////////////////////////
1273/// Determine total error matrix on the vector Ax.
1274
1276{
1278
1279 // errors from input vector and from background subtraction
1280 TMatrixDSparse *emat_sum=new TMatrixDSparse(*fVyy);
1281
1282 // uncorrelated systematic error
1283 if(fEmatUncorrAx) {
1284 AddMSparse(emat_sum,1.0,fEmatUncorrAx);
1285 }
1286 TMapIter sysErrPtr(fDeltaCorrAx);
1287 const TObject *key;
1288
1289 // correlated systematic errors
1290 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1291 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1293 AddMSparse(emat_sum,1.0,emat);
1294 DeleteMatrix(&emat);
1295 }
1296 // error on tau
1297 if(fDeltaSysTau) {
1299 TMatrixDSparse *emat_tau=
1300 MultiplyMSparseMSparseTranspVector(Adx_tau,Adx_tau,0);
1301 DeleteMatrix(&Adx_tau);
1302 AddMSparse(emat_sum,1.0,emat_tau);
1303 DeleteMatrix(&emat_tau);
1304 }
1305 return emat_sum;
1306}
1307
1308////////////////////////////////////////////////////////////////////////////////
1309/// Determine total error matrix on the vector x.
1310
1312{
1314
1315 // errors from input vector and from background subtraction
1316 TMatrixDSparse *emat_sum=new TMatrixDSparse(*GetVxx());
1317
1318 // uncorrelated systematic error
1319 if(fEmatUncorrX) {
1320 AddMSparse(emat_sum,1.0,fEmatUncorrX);
1321 }
1322 TMapIter sysErrPtr(fDeltaCorrX);
1323 const TObject *key;
1324
1325 // correlated systematic errors
1326 for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1327 const TMatrixDSparse *delta=(TMatrixDSparse *)((const TPair *)*sysErrPtr)->Value();
1329 AddMSparse(emat_sum,1.0,emat);
1330 DeleteMatrix(&emat);
1331 }
1332 // error on tau
1333 if(fDeltaSysTau) {
1334 TMatrixDSparse *emat_tau=
1336 AddMSparse(emat_sum,1.0,emat_tau);
1337 DeleteMatrix(&emat_tau);
1338 }
1339 return emat_sum;
1340}
1341
1342
1343////////////////////////////////////////////////////////////////////////////////
1344/// Calculate total chi**2 including all systematic errors.
1345
1347{
1348
1350
1351 Int_t rank=0;
1352 TMatrixDSparse *v=InvertMSparseSymmPos(emat_sum,&rank);
1355 DeleteMatrix(&v);
1356 Double_t r=0.0;
1357 const Int_t *vdy_rows=vdy->GetRowIndexArray();
1358 const Double_t *vdy_data=vdy->GetMatrixArray();
1359 for(Int_t i=0;i<vdy->GetNrows();i++) {
1360 if(vdy_rows[i+1]>vdy_rows[i]) {
1361 r += vdy_data[vdy_rows[i]] * dy(i,0);
1362 }
1363 }
1364 DeleteMatrix(&vdy);
1365 DeleteMatrix(&emat_sum);
1366 return r;
1367}
1368
1369////////////////////////////////////////////////////////////////////////////////
1370/// Get global correlatiocn coefficients, summing up all contributions.
1371///
1372/// \param[out] rhoi histogram which will be filled
1373/// \param[in] binMap (default=0) remapping of histogram bins
1374/// \param[out] invEmat (default=0) inverse of error matrix
1375///
1376/// return the global correlation coefficients, including all error
1377/// sources. If <b>invEmat</b> is nonzero, the inverse of the error
1378/// matrix is returned in that histogram
1379///
1380/// the array <b>binMap</b> is explained with the method GetOutput().
1381///
1382/// get global correlation coefficients including systematic,statistical,background,tau errors
1383/// - rhoi: output histogram
1384/// - binMap: for each global bin, indicate in which histogram bin
1385/// to store its content
1386/// - invEmat: output histogram for inverse of error matrix
1387/// (pointer may zero if inverse is not requested)
1388
1389void TUnfoldSys::GetRhoItotal(TH1 *rhoi,const Int_t *binMap,TH2 *invEmat)
1390{
1391 ClearHistogram(rhoi,-1.);
1393 GetRhoIFromMatrix(rhoi,emat_sum,binMap,invEmat);
1394
1395 DeleteMatrix(&emat_sum);
1396}
1397
1398////////////////////////////////////////////////////////////////////////////////
1399/// Scale columns of a matrix by the corresponding rows of a vector.
1400///
1401/// \param[inout] m matrix
1402/// \param[in] v vector
1403///
1404/// the entries m<sub>ij</sub> are multiplied by v<sub>j</sub>.
1405///
1406/// scale columns of m by the corresponding rows of v
1407/// input:
1408/// - m: pointer to sparse matrix of dimension NxM
1409/// - v: pointer to matrix of dimension Mx1
1410
1413{
1414 if((m->GetNcols() != v->GetNrows())||(v->GetNcols()!=1)) {
1415 Fatal("ScaleColumnsByVector error",
1416 "matrix cols/vector rows %d!=%d OR vector cols %d !=1\n",
1417 m->GetNcols(),v->GetNrows(),v->GetNcols());
1418 }
1419 const Int_t *rows_m=m->GetRowIndexArray();
1420 const Int_t *cols_m=m->GetColIndexArray();
1421 Double_t *data_m=m->GetMatrixArray();
1422 const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
1423 if(v_sparse) {
1424 const Int_t *rows_v=v_sparse->GetRowIndexArray();
1425 const Double_t *data_v=v_sparse->GetMatrixArray();
1426 for(Int_t i=0;i<m->GetNrows();i++) {
1427 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1428 Int_t j=cols_m[index_m];
1429 Int_t index_v=rows_v[j];
1430 if(index_v<rows_v[j+1]) {
1431 data_m[index_m] *= data_v[index_v];
1432 } else {
1433 data_m[index_m] =0.0;
1434 }
1435 }
1436 }
1437 } else {
1438 for(Int_t i=0;i<m->GetNrows();i++) {
1439 for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1440 Int_t j=cols_m[index_m];
1441 data_m[index_m] *= (*v)(j,0);
1442 }
1443 }
1444 }
1445}
1446
1447////////////////////////////////////////////////////////////////////////////////
1448/// Map delta to hist_delta, possibly summing up bins.
1449///
1450/// \param[out] hist_delta result histogram
1451/// \param[in] delta vector to be mapped to the histogram
1452/// \param[in] binMap mapping of histogram bins
1453///
1454/// groups of bins of <b>delta</b> are mapped to bins of
1455/// <b>hist_delta</b>. The histogram contents are set to the sum over
1456/// the group of bins. The histogram errors are reset to zero.
1457///
1458/// The array <b>binMap</b> is explained with the method GetOutput()
1459///
1460/// sum over bins of *delta, as defined in binMap,fXToHist
1461/// - hist_delta: histogram to return summed vector
1462/// - delta: vector to sum and remap
1463
1465(TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap)
1466{
1467 Int_t nbin=hist_delta->GetNbinsX();
1468 Double_t *c=new Double_t[nbin+2];
1469 for(Int_t i=0;i<nbin+2;i++) {
1470 c[i]=0.0;
1471 }
1472 if(delta) {
1473 Int_t binMapSize = fHistToX.GetSize();
1474 const Double_t *delta_data=delta->GetMatrixArray();
1475 const Int_t *delta_rows=delta->GetRowIndexArray();
1476 for(Int_t i=0;i<binMapSize;i++) {
1477 Int_t destBinI=binMap ? binMap[i] : i;
1478 Int_t srcBinI=fHistToX[i];
1479 if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
1480 Int_t index=delta_rows[srcBinI];
1481 if(index<delta_rows[srcBinI+1]) {
1482 c[destBinI]+=delta_data[index];
1483 }
1484 }
1485 }
1486 }
1487 for(Int_t i=0;i<nbin+2;i++) {
1488 hist_delta->SetBinContent(i,c[i]);
1489 hist_delta->SetBinError(i,0.0);
1490 }
1491 delete[] c;
1492}
1493
1494////////////////////////////////////////////////////////////////////////////////
1495/// Get a new list of all systematic uuncertainty sources.
1496///
1497/// The user is responsible for deleting the list
1498/// get list of names of systematic sources
1499
1501 TSortedList *r=new TSortedList();
1502 TMapIter i(fSysIn);
1503 for(const TObject *key=i.Next();key;key=i.Next()) {
1504 r->Add(((TObjString *)key)->Clone());
1505 }
1506 return r;
1507}
1508
1509////////////////////////////////////////////////////////////////////////////////
1510/// Get a new list of all background sources.
1511///
1512/// The user is responsible for deleting the list
1513/// get list of name of background sources
1514
1516 TSortedList *r=new TSortedList();
1517 TMapIter i(fBgrIn);
1518 for(const TObject *key=i.Next();key;key=i.Next()) {
1519 r->Add(((TObjString *)key)->Clone());
1520 }
1521 return r;
1522}
ROOT::R::TRInterface & r
Definition: Object.C:4
#define c(i)
Definition: RSha256.hxx:101
const Bool_t kFALSE
Definition: RtypesCore.h:90
double Double_t
Definition: RtypesCore.h:57
const Bool_t kTRUE
Definition: RtypesCore.h:89
#define ClassImp(name)
Definition: Rtypes.h:361
char name[80]
Definition: TGX11.cxx:109
TMatrixTSparse< Double_t > TMatrixDSparse
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:22
Int_t GetSize() const
Definition: TArray.h:47
virtual Int_t GetEntries() const
Definition: TCollection.h:177
The TH1 histogram class.
Definition: TH1.h:56
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8519
virtual Int_t GetNbinsX() const
Definition: TH1.h:292
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:8662
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:8678
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4907
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:88
Iterator of map.
Definition: TMap.h:147
TObject * Next()
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)
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
Check if a (key,value) pair exists with keyname as name of the key.
Definition: TMap.cxx:215
void Clear(Option_t *option="")
Remove all (key,value) pairs from the map.
Definition: TMap.cxx:97
Int_t GetNrows() const
Definition: TMatrixTBase.h:124
virtual const Int_t * GetRowIndexArray() const
virtual const Element * GetMatrixArray() const
virtual const Int_t * GetColIndexArray() const
Collectable string class.
Definition: TObjString.h:28
const TString & GetString() const
Definition: TObjString.h:46
Mother of all ROOT objects.
Definition: TObject.h:37
virtual TObject * Clone(const char *newname="") const
Make a clone of an object using the Streamer facility.
Definition: TObject.cxx:144
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:877
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:891
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:919
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:865
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:131
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:418
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
void GetEmatrixSysSource(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance contribution from a systematic variation of the response matrix.
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.
Double_t GetChi2Sys(void)
Calculate total chi**2 including all systematic errors.
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.
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 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=0, const Int_t *binMap=0, 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
void GetEmatrixSysBackgroundScale(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance contribution from background normalisation uncertainty.
TMatrixDSparse * fVyyData
Input: error on fY prior to bgr subtraction.
Definition: TUnfoldSys.h:78
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=0)
Correlated one-sigma shifts correspinding to a given systematic uncertainty.
Definition: TUnfoldSys.cxx:999
TMatrixDSparse * fEmatUncorrAx
Result: syst.error from fDA2 on fAx.
Definition: TUnfoldSys.h:82
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
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance matrix contribution from input measurement uncertainties.
Bool_t GetDeltaSysBackgroundScale(TH1 *delta, const char *source, const Int_t *binMap=0)
Correlated one-sigma shifts from background normalisation uncertainty.
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
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=0)
Correlated one-sigma shifts from shifting tau.
virtual TMatrixDSparse * PrepareUncorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2)
Propagate uncorrelated systematic errors to a covariance matrix.
Definition: TUnfoldSys.cxx:838
virtual void ClearResults(void)
Clear all data members which depend on the unfolding results.
Definition: TUnfoldSys.cxx:616
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=0, const TH2 *hist_vyy_inv=0)
Define the input data for subsequent calls to DoUnfold(Double_t).
Definition: TUnfoldSys.cxx:444
TMatrixDSparse * fEmatUncorrX
Result: syst.error from fDA2 on fX.
Definition: TUnfoldSys.h:80
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 GetRhoItotal(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0)
Get global correlatiocn coefficients, summing up all contributions.
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 GetEmatrixSysTau(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance matrix contribution from error on regularisation parameter.
void AddSysError(const TH2 *sysError, const char *name, EHistMap histmap, ESysErrMode mode)
Specify a correlated systematic uncertainty.
Definition: TUnfoldSys.cxx:248
virtual ~TUnfoldSys(void)
Definition: TUnfoldSys.cxx:119
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
void GetEmatrixSysUncorr(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance contribution from uncorrelated uncertainties of the response matrix.
Definition: TUnfoldSys.cxx:731
void GetEmatrixTotal(TH2 *ematrix, const Int_t *binMap=0)
Get total error matrix, summing up all contributions.
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance contribution from background uncorrelated uncertainty.
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:773
TMatrixDSparse * MultiplyMSparseTranspMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply a transposed Sparse matrix with another sparse matrix,.
Definition: TUnfold.cxx:693
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:833
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:592
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:195
void ClearHistogram(TH1 *h, Double_t x=0.) const
Initialize bin contents and bin errors for a given histogram.
Definition: TUnfold.cxx:3643
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
virtual Int_t SetInput(const TH1 *hist_y, Double_t scaleBias=0.0, Double_t oneOverZeroError=0.0, const TH2 *hist_vyy=0, const TH2 *hist_vyy_inv=0)
Define input data for subsequent calls to DoUnfold(tau).
Definition: TUnfold.cxx:2300
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=0) const
Get output covariance matrix, possibly cumulated over several bins.
Definition: TUnfold.cxx:3410
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:3343
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:1008
Double_t fTauSquared
regularisation parameter tau squared
Definition: TUnfold.h:162
virtual void ClearResults(void)
Reset all results.
Definition: TUnfold.cxx:218
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Multiply two sparse matrices.
Definition: TUnfold.cxx:617
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:930
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:3528
const TMatrixDSparse * GetDXDY(void) const
matrix of derivatives dx/dy
Definition: TUnfold.h:246
TMatrixDSparse * fA
response matrix A
Definition: TUnfold.h:150
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)
Definition: TMath.h:681
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
const char * Value
Definition: TXMLSetup.cxx:72
auto * m
Definition: textangle.C:8
static long int sum(long int i)
Definition: Factory.cxx:2275