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