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