ROOT  6.06/09
Reference Guide
TUnfoldSys.cxx
Go to the documentation of this file.
1 // Author: Stefan Schmitt
2 // DESY, 23/01/09
3 
4 // Version 17.1, bug fix with background uncertainty
5 //
6 // History:
7 // Version 17.0, possibility to specify an error matrix with SetInput
8 // Version 16.1, parallel to changes in TUnfold
9 // Version 16.0, parallel to changes in TUnfold
10 // Version 15, fix bugs with uncorr. uncertainties, add backgnd subtraction
11 // Version 14, remove some print-out, do not add unused sys.errors
12 // Version 13, support for systematic errors
13 
14 /** \class TUnfoldSys
15  \ingroup Hist
16  TUnfold is used to decompose a measurement y into several sources x
17  given the measurement uncertainties and a matrix of migrations A
18 
19  TUnfoldSys adds error propagation of systematic errors to TUnfold
20  Also, background sources (with errors) can be subtracted.
21 
22  **For most applications, it is better to use TUnfoldDensity
23  instead of using TUnfoldSys or TUnfold**
24 
25  If you use this software, please consider the following citation
26  S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]
27 
28  More documentation and updates are available on
29  http://www.desy.de/~sschmitt
30 
31  The following sources of systematic error are considered in TUnfoldSys
32 
33  (a) uncorrelated errors on the input matrix histA, taken as the
34  errors provided with the histogram.
35  These are typically statistical errors from finite Monte Carlo samples
36 
37  (b) correlated shifts of the input matrix histA. These shifts are taken
38  as one-sigma effects when switchig on a given error soure.
39  several such error sources may be defined
40 
41  (c) a systematic error on the regularisation parameter tau
42 
43  (d) uncorrelated errors on background sources, taken as the errors
44  provided with the background histograms
45 
46  (e) scale errors on background sources
47 
48  In addition there is the (statistical) uncertainty of the input vector (i)
49 
50  Source (a) is providede with the original histogram histA
51  TUnfoldSys(histA,...)
52 
53  Sources (b) are added by calls to
54  AddSysError()
55 
56  The systematic uncertainty on tau (c) is set by
57  SetTauError()
58 
59  Backgound sources causing errors of type (d) and (e) are added by
60  SubtractBackground()
61 
62  NOTE:
63  Systematic errors (a), (b), (c) are propagated to the result
64  AFTER unfolding
65 
66  Background errors (d) and (e) are added to the data errors
67  BEFORE unfolding
68 
69  For this reason:
70  errors of type (d) and (e) are INCLUDED in the standard error matrix
71  and other methods provided by the base class TUnfold:
72  GetOutput()
73  GetEmatrix()
74  ...
75  whereas errors of type (a), (b), (c) are NOT INCLUDED in the methods
76  provided by the base class TUnfold.
77 
78  ## Accessing error matrices:
79  The error sources (b),(c) and (e) propagate to shifts of the result.
80  These shifts may be accessed as histograms using the methods
81  GetDeltaSysSource() corresponds to (b)
82  GetDeltaSysTau() corresponds to (c)
83  GetDeltaSysBackgroundScale() corresponds to (e)
84  The error sources (a) and (d) originate from many uncorrelated errors,
85  which in general are NOT uncorrelated on the result vector.
86  Thus, there is no corresponding shift of the output vector, only error
87  matrices are available
88 
89  Method to get error matrix | corresponds to error sources
90  -----------------------------|---------------------------------
91  GetEmatrixSysUncorr() | (a)
92  GetEmatrixSysSource() | (b)
93  GetEmatrixSysTau() | (c)
94  GetEmatrixSysBackgroundUncorr() | (d)
95  GetEmatrixSysBackgroundScale() | (e)
96  GetEmatrixInput() | (i)
97  GetEmatrix() | (i)+(d)+(e)
98  GetEmatrixTotal() | (i)+(a)+(b)+(c)+(d)+(e)
99 
100  Error matrices can be added to existing histograms.
101  This is useful to retreive the sum of several error matrices.
102  If the last argument of the GetEmatrixXXX() methods is set to kFALSE,
103  the histogram is not cleared, but the error matrix is simply added to the
104  existing histogram
105 */
106 
107 /*
108  This file is part of TUnfold.
109 
110  TUnfold is free software: you can redistribute it and/or modify
111  it under the terms of the GNU General Public License as published by
112  the Free Software Foundation, either version 3 of the License, or
113  (at your option) any later version.
114 
115  TUnfold is distributed in the hope that it will be useful,
116  but WITHOUT ANY WARRANTY; without even the implied warranty of
117  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
118  GNU General Public License for more details.
119 
120  You should have received a copy of the GNU General Public License
121  along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
122 */
123 
124 #include <iostream>
125 #include <TMap.h>
126 #include <TMath.h>
127 #include <TObjString.h>
128 #include <RVersion.h>
129 #include <cmath>
130 
131 #include "TUnfoldSys.h"
132 
134 
135 TUnfoldSys::TUnfoldSys(void)
136 {
137  // set all pointers to zero
138  InitTUnfoldSys();
139 }
140 
142 (const TH2 *hist_A, EHistMap histmap, ERegMode regmode,EConstraint constraint)
143  : TUnfold(hist_A,histmap,regmode,constraint)
144 {
145  // arguments:
146  // hist_A: matrix that describes the migrations
147  // histmap: mapping of the histogram axes to the unfolding output
148  // regmode: global regularisation mode
149  // data members initialized to something different from zero:
150  // fDA2, fDAcol
151 
152  // initialize TUnfold
153  InitTUnfoldSys();
154 
155  // svae underflow and overflow bins
156  fAoutside = new TMatrixD(GetNx(),2);
157  // save the romalized errors on hist_A
158  // to the matrices fDAinRelSq and fDAinColRelSq
159  fDAinColRelSq = new TMatrixD(GetNx(),1);
160 
161  Int_t nmax=GetNx()*GetNy();
162  Int_t *rowDAinRelSq = new Int_t[nmax];
163  Int_t *colDAinRelSq = new Int_t[nmax];
164  Double_t *dataDAinRelSq = new Double_t[nmax];
165 
166  Int_t da_nonzero=0;
167  for (Int_t ix = 0; ix < GetNx(); ix++) {
168  Int_t ibinx = fXToHist[ix];
169  Double_t sum_sq= fSumOverY[ix]*fSumOverY[ix];
170  for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
171  Double_t dz;
172  if (histmap == kHistMapOutputHoriz) {
173  dz = hist_A->GetBinError(ibinx, ibiny);
174  } else {
175  dz = hist_A->GetBinError(ibiny, ibinx);
176  }
177  Double_t normerr_sq=dz*dz/sum_sq;
178  // quadratic sum of all errors from all bins,
179  // including under/overflow bins
180  (*fDAinColRelSq)(ix,0) += normerr_sq;
181 
182  if(ibiny==0) {
183  // underflow bin
184  if (histmap == kHistMapOutputHoriz) {
185  (*fAoutside)(ix,0)=hist_A->GetBinContent(ibinx, ibiny);
186  } else {
187  (*fAoutside)(ix,0)=hist_A->GetBinContent(ibiny, ibinx);
188  }
189  } else if(ibiny==GetNy()+1) {
190  // overflow bins
191  if (histmap == kHistMapOutputHoriz) {
192  (*fAoutside)(ix,1)=hist_A->GetBinContent(ibinx, ibiny);
193  } else {
194  (*fAoutside)(ix,1)=hist_A->GetBinContent(ibiny, ibinx);
195  }
196  } else {
197  // error on this bin
198  rowDAinRelSq[da_nonzero]=ibiny-1;
199  colDAinRelSq[da_nonzero] = ix;
200  dataDAinRelSq[da_nonzero] = normerr_sq;
201  if(dataDAinRelSq[da_nonzero]>0.0) da_nonzero++;
202  }
203  }
204  }
205  if(da_nonzero) {
206  fDAinRelSq = CreateSparseMatrix(GetNy(),GetNx(),da_nonzero,
207  rowDAinRelSq,colDAinRelSq,dataDAinRelSq);
208  } else {
209  DeleteMatrix(&fDAinColRelSq);
210  }
211  delete[] rowDAinRelSq;
212  delete[] colDAinRelSq;
213  delete[] dataDAinRelSq;
214 }
215 
217 (const TH2 *sysError,const char *name,EHistMap histmap,ESysErrMode mode)
218 {
219  // add a correlated error source
220  // sysError: alternative matrix or matrix of absolute/relative shifts
221  // name: name of the error source
222  // histmap: mapping of the histogram axes to the unfolding output
223  // mode: format of the error source
224 
225  if(fSysIn->FindObject(name)) {
226  Error("AddSysError","Source %s given twice, ignoring 2nd call.\n",name);
227  } else {
228  // a copy of fA is made. It can be accessed inside the loop
229  // without having to take care that the sparse structure of *fA
230  // may be accidentally destroyed by asking for an element which is zero.
231  TMatrixD aCopy(*fA);
232 
233  Int_t nmax= GetNx()*GetNy();
234  Double_t *data=new Double_t[nmax];
235  Int_t *cols=new Int_t[nmax];
236  Int_t *rows=new Int_t[nmax];
237  nmax=0;
238  for (Int_t ix = 0; ix < GetNx(); ix++) {
239  Int_t ibinx = fXToHist[ix];
240  Double_t sum=0.0;
241  for(Int_t loop=0;loop<2;loop++) {
242  for (Int_t ibiny = 0; ibiny <= GetNy()+1; ibiny++) {
243  Double_t z;
244  // get bin content, depending on histmap
245  if (histmap == kHistMapOutputHoriz) {
246  z = sysError->GetBinContent(ibinx, ibiny);
247  } else {
248  z = sysError->GetBinContent(ibiny, ibinx);
249  }
250  // correct to absolute numbers
251  if(mode!=kSysErrModeMatrix) {
252  Double_t z0;
253  if((ibiny>0)&&(ibiny<=GetNy())) {
254  z0=aCopy(ibiny-1,ix)*fSumOverY[ix];
255  } else if(ibiny==0) {
256  z0=(*fAoutside)(ix,0);
257  } else {
258  z0=(*fAoutside)(ix,1);
259  }
260  if(mode==kSysErrModeShift) {
261  z += z0;
262  } else if(mode==kSysErrModeRelative) {
263  z = z0*(1.+z);
264  }
265  }
266  if(loop==0) {
267  // sum up all entries, including overflow bins
268  sum += z;
269  } else {
270  if((ibiny>0)&&(ibiny<=GetNy())) {
271  // save normalized matrix of shifts,
272  // excluding overflow bins
273  rows[nmax]=ibiny-1;
274  cols[nmax]=ix;
275  if(sum>0.0) {
276  data[nmax]=z/sum-aCopy(ibiny-1,ix);
277  } else {
278  data[nmax]=0.0;
279  }
280  if(data[nmax] != 0.0) nmax++;
281  }
282  }
283  }
284  }
285  }
286  if(nmax==0) {
287  Error("AddSysError",
288  "source %s has no influence and has not been added.\n",name);
289  } else {
290  TMatrixDSparse *dsys=CreateSparseMatrix(GetNy(),GetNx(),
291  nmax,rows,cols,data);
292  fSysIn->Add(new TObjString(name),dsys);
293  }
294  delete[] data;
295  delete[] rows;
296  delete[] cols;
297  }
298 }
299 
301 {
302  // performs background subtraction
303  // fY = fYData - fBgrIn
304  // fVyy = fVyyData + fBgrErrUncorr^2 + fBgrErrCorr * fBgrErrCorr#
305  // fVyyinv = fVyy^(-1)
306 
307  if(fYData) {
308  DeleteMatrix(&fY);
309  DeleteMatrix(&fVyy);
310  if(fBgrIn->GetEntries()<=0) {
311  // simple copy
312  fY=new TMatrixD(*fYData);
314  } else {
315  if(GetVyyInv()) {
316  Warning("DoBackgroundSubtraction",
317  "inverse error matrix from user input,"
318  " not corrected for background");
319  }
320  // copy of the data
321  fY=new TMatrixD(*fYData);
322  // subtract background from fY
323  const TObject *key;
324  {
325  TMapIter bgrPtr(fBgrIn);
326  for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
327  const TMatrixD *bgr=(const TMatrixD *)
328 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
329  ((const TPair *)*bgrPtr)->Value()
330 #else
331  fBgrIn->GetValue(((const TObjString *)key)->GetString())
332 #endif
333  ;
334  for(Int_t i=0;i<GetNy();i++) {
335  (*fY)(i,0) -= (*bgr)(i,0);
336  }
337  }
338  }
339  // copy original error matrix
340  TMatrixD vyy(*fVyyData);
341  // determine used bins
343  const Int_t *vyydata_rows=fVyyData->GetRowIndexArray();
344  const Int_t *vyydata_cols=fVyyData->GetColIndexArray();
345  const Double_t *vyydata_data=fVyyData->GetMatrixArray();
346  Int_t *usedBin=new Int_t[ny];
347  for(Int_t i=0;i<ny;i++) {
348  usedBin[i]=0;
349  }
350  for(Int_t i=0;i<ny;i++) {
351  for(Int_t k=vyydata_rows[i];k<vyydata_rows[i+1];k++) {
352  if(vyydata_data[k]>0.0) {
353  usedBin[i]++;
354  usedBin[vyydata_cols[k]]++;
355  }
356  }
357  }
358  // add uncorrelated background errors
359  {
360  TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
361  for(key=bgrErrUncorrSqPtr.Next();key;
362  key=bgrErrUncorrSqPtr.Next()) {
363  const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)
364 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
365  ((const TPair *)*bgrErrUncorrSqPtr)->Value()
366 #else
367  fBgrErrUncorrInSq->GetValue(((const TObjString *)key)
368  ->GetString())
369 #endif
370  ;
371  for(Int_t yi=0;yi<ny;yi++) {
372  if(!usedBin[yi]) continue;
373  vyy(yi,yi) +=(*bgrerruncorrSquared)(yi,0);
374  }
375  }
376  }
377  // add correlated background errors
378  {
379  TMapIter bgrErrScalePtr(fBgrErrScaleIn);
380  for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
381  const TMatrixD *bgrerrscale=(const TMatrixD *)
382 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
383  ((const TPair *)*bgrErrScalePtr)->Value()
384 #else
385  fBgrErrScaleIn->GetValue(((const TObjString *)key)
386  ->GetString())
387 #endif
388  ;
389  for(Int_t yi=0;yi<ny;yi++) {
390  if(!usedBin[yi]) continue;
391  for(Int_t yj=0;yj<ny;yj++) {
392  if(!usedBin[yj]) continue;
393  vyy(yi,yj) +=(*bgrerrscale)(yi,0)* (*bgrerrscale)(yj,0);
394  }
395  }
396  }
397  }
398  delete[] usedBin;
399  usedBin=0;
400 
401  // convert to sparse matrix
402  fVyy=new TMatrixDSparse(vyy);
403 
404  }
405  } else {
406  Fatal("DoBackgroundSubtraction","No input vector defined");
407  }
408 }
409 
410 Int_t TUnfoldSys::SetInput(const TH1 *hist_y,Double_t scaleBias,
411  Double_t oneOverZeroError,const TH2 *hist_vyy,
412  const TH2 *hist_vyy_inv)
413 {
414  // Define the input data for subsequent calls to DoUnfold(Double_t)
415  // input: input distribution with errors
416  // scaleBias: scale factor applied to the bias
417  // oneOverZeroError: for bins with zero error, this number defines 1/error.
418  // Return value: number of bins with bad error
419  // +10000*number of unconstrained output bins
420  // Note: return values>=10000 are fatal errors,
421  // for the given input, the unfolding can not be done!
422  // Calls the SetInput method of the base class, then renames the input
423  // vectors fY and fVyy, then performs the background subtraction
424  // Data members modified:
425  // fYData,fY,fVyyData,fVyy,fVyyinvData,fVyyinv
426  // and those modified by TUnfold::SetInput()
427  // and those modified by DoBackgroundSubtraction()
428 
429  Int_t r=TUnfold::SetInput(hist_y,scaleBias,oneOverZeroError,hist_vyy,
430  hist_vyy_inv);
431  fYData=fY;
432  fY=0;
433  fVyyData=fVyy;
434  fVyy=0;
436 
437  return r;
438 }
439 
441 (const TH1 *bgr,const char *name,Double_t scale,Double_t scale_error)
442 {
443  // Store background source
444  // bgr: background distribution with uncorrelated errors
445  // name: name of this background source
446  // scale: scale factor applied to the background
447  // scaleError: error on scale factor (correlated error)
448  //
449  // Data members modified:
450  // fBgrIn,fBgrErrUncorrInSq,fBgrErrScaleIn
451  // and those modified by DoBackgroundSubtraction()
452 
453  // save background source
454  if(fBgrIn->FindObject(name)) {
455  Error("SubtractBackground","Source %s given twice, ignoring 2nd call.\n",
456  name);
457  } else {
458  TMatrixD *bgrScaled=new TMatrixD(GetNy(),1);
459  TMatrixD *bgrErrUncSq=new TMatrixD(GetNy(),1);
460  TMatrixD *bgrErrCorr=new TMatrixD(GetNy(),1);
461  for(Int_t row=0;row<GetNy();row++) {
462  (*bgrScaled)(row,0) = scale*bgr->GetBinContent(row+1);
463  (*bgrErrUncSq)(row,0) =
464  TMath::Power(scale*bgr->GetBinError(row+1),2.);
465  (*bgrErrCorr)(row,0) = scale_error*bgr->GetBinContent(row+1);
466  }
467  fBgrIn->Add(new TObjString(name),bgrScaled);
468  fBgrErrUncorrInSq->Add(new TObjString(name),bgrErrUncSq);
469  fBgrErrScaleIn->Add(new TObjString(name),bgrErrCorr);
470  if(fYData) {
471  DoBackgroundSubtraction();
472  } else {
473  Info("SubtractBackground",
474  "Background subtraction prior to setting input data");
475  }
476  }
477 }
478 
480 (TH1 *bgrHist,const char *bgrSource,const Int_t *binMap,
481  Int_t includeError,Bool_t clearHist) const
482 {
483  if(clearHist) {
484  ClearHistogram(bgrHist);
485  }
486  // add all background sources
487  const TObject *key;
488  {
489  TMapIter bgrPtr(fBgrIn);
490  for(key=bgrPtr.Next();key;key=bgrPtr.Next()) {
491  TString bgrName=((const TObjString *)key)->GetString();
492  if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
493  const TMatrixD *bgr=(const TMatrixD *)
494 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
495  ((const TPair *)*bgrPtr)->Value()
496 #else
497  fBgrIn->GetValue(bgrName)
498 #endif
499  ;
500  for(Int_t i=0;i<GetNy();i++) {
501  Int_t destBin=binMap[i];
502  bgrHist->SetBinContent(destBin,bgrHist->GetBinContent(destBin)+
503  (*bgr)(i,0));
504  }
505  }
506  }
507  // add uncorrelated background errors
508  if(includeError &1) {
509  TMapIter bgrErrUncorrSqPtr(fBgrErrUncorrInSq);
510  for(key=bgrErrUncorrSqPtr.Next();key;key=bgrErrUncorrSqPtr.Next()) {
511  TString bgrName=((const TObjString *)key)->GetString();
512  if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
513  const TMatrixD *bgrerruncorrSquared=(TMatrixD const *)
514 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
515  ((const TPair *)*bgrErrUncorrSqPtr)->Value()
516 #else
517  fBgrErrUncorrInSq->GetValue(((const TObjString *)key)
518  ->GetString())
519 #endif
520  ;
521  for(Int_t i=0;i<GetNy();i++) {
522  Int_t destBin=binMap[i];
523  bgrHist->SetBinError
524  (destBin,TMath::Sqrt
525  ((*bgrerruncorrSquared)(i,0)+
526  TMath::Power(bgrHist->GetBinError(destBin),2.)));
527  }
528  }
529  }
530  if(includeError & 2) {
531  TMapIter bgrErrScalePtr(fBgrErrScaleIn);
532  for(key=bgrErrScalePtr.Next();key;key=bgrErrScalePtr.Next()) {
533  TString bgrName=((const TObjString *)key)->GetString();
534  if(bgrSource && bgrName.CompareTo(bgrSource)) continue;
535  const TMatrixD *bgrerrscale=(TMatrixD const *)
536 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
537  ((const TPair *)*bgrErrScalePtr)->Value()
538 #else
539  fBgrErrScaleIn->GetValue(((const TObjString *)key)->GetString())
540 #endif
541  ;
542  for(Int_t i=0;i<GetNy();i++) {
543  Int_t destBin=binMap[i];
544  bgrHist->SetBinError(destBin,hypot((*bgrerrscale)(i,0),
545  bgrHist->GetBinError(destBin)));
546  }
547  }
548  }
549 }
550 
552 {
553  // initialize pointers and TMaps
554 
555  // input
556  fDAinRelSq = 0;
557  fDAinColRelSq = 0;
558  fAoutside = 0;
559  fBgrIn = new TMap();
560  fBgrErrUncorrInSq = new TMap();
561  fBgrErrScaleIn = new TMap();
562  fSysIn = new TMap();
563 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
568 #else
569  fBgrIn->SetOwner();
572  fSysIn->SetOwner();
573 #endif
574  // results
575  fEmatUncorrX = 0;
576  fEmatUncorrAx = 0;
577  fDeltaCorrX = new TMap();
578  fDeltaCorrAx = new TMap();
579 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
582 #else
585 #endif
586  fDeltaSysTau = 0;
587  fDtau=0.0;
588  fYData=0;
589  fVyyData=0;
590 }
591 
593 {
594  // delete all data members
597  delete fBgrIn;
598  delete fBgrErrUncorrInSq;
599  delete fBgrErrScaleIn;
600  delete fSysIn;
601  ClearResults();
602  delete fDeltaCorrX;
603  delete fDeltaCorrAx;
606 }
607 
609 {
610  // clear all data members which depend on the unfolding results
614  fDeltaCorrX->Clear();
615  fDeltaCorrAx->Clear();
617 }
618 
620 {
621  // calculations required for syst.error
622  // data members modified
623  // fEmatUncorrX, fEmatUncorrAx, fDeltaCorrX, fDeltaCorrAx
624  if(!fEmatUncorrX) {
626  }
627  TMatrixDSparse *AM0=0,*AM1=0;
628  if(!fEmatUncorrAx) {
629  if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
630  if(!AM1) {
632  Int_t *rows_cols=new Int_t[GetNy()];
633  Double_t *data=new Double_t[GetNy()];
634  for(Int_t i=0;i<GetNy();i++) {
635  rows_cols[i]=i;
636  data[i]=1.0;
637  }
639  (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
640  delete[] data;
641  delete[] rows_cols;
642  AddMSparse(AM1,-1.,one);
643  DeleteMatrix(&one);
645  }
646  }
647  if((!fDeltaSysTau )&&(fDtau>0.0)) {
652  for(Int_t i=0;i<n;i++) {
653  data[i] *= scale;
654  }
655  }
656 
657  TMapIter sysErrIn(fSysIn);
658  const TObjString *key;
659 
660  // calculate individual systematic errors
661  for(key=(const TObjString *)sysErrIn.Next();key;
662  key=(const TObjString *)sysErrIn.Next()) {
663 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
664  const TMatrixDSparse *dsys=
665  (const TMatrixDSparse *)((const TPair *)*sysErrIn)->Value();
666 #else
667  const TMatrixDSparse *dsys=
668  (const TMatrixDSparse *)(fSysIn->GetValue(key->GetString()));
669 #endif
670  const TPair *named_emat=(const TPair *)
672  if(!named_emat) {
674  fDeltaCorrX->Add(new TObjString(*key),emat);
675  }
676  named_emat=(const TPair *)fDeltaCorrAx->FindObject(key->GetString());
677  if(!named_emat) {
678  if(!AM0) AM0=MultiplyMSparseMSparse(fA,GetDXDAM(0));
679  if(!AM1) {
681  Int_t *rows_cols=new Int_t[GetNy()];
682  Double_t *data=new Double_t[GetNy()];
683  for(Int_t i=0;i<GetNy();i++) {
684  rows_cols[i]=i;
685  data[i]=1.0;
686  }
688  (GetNy(),GetNy(),GetNy(),rows_cols, rows_cols,data);
689  delete[] data;
690  delete[] rows_cols;
691  AddMSparse(AM1,-1.,one);
692  DeleteMatrix(&one);
694  }
695  TMatrixDSparse *emat=PrepareCorrEmat(AM0,AM1,dsys);
696  fDeltaCorrAx->Add(new TObjString(*key),emat);
697  }
698  }
699  DeleteMatrix(&AM0);
700  DeleteMatrix(&AM1);
701 }
702 
704 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
705 {
706  // get output error contribution from statistical fluctuations in A
707  // ematrix: output error matrix histogram
708  // binMap: see method GetEmatrix()
709  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
710  // data members modified:
711  // fVYAx, fESparse, fEAtV, fErrorAStat
712  PrepareSysError();
713  ErrorMatrixToHist(ematrix,fEmatUncorrX,binMap,clearEmat);
714 }
715 
717 (const TMatrixDSparse *m_0,const TMatrixDSparse *m_1)
718 {
719  // propagate uncorrelated systematic errors to a covariance matrix
720  // m0,m1 : coefficients (matrices) for propagating the errors
721  //
722  // the error matrix is calculated by standard error propagation, where the
723  // derivative of the result vector X wrt the matrix A is given by
724  //
725  // dX_k / dA_ij = M0_kj * Z0_i - M1_ki * Z1_j
726  //
727  // where:
728  // the matrices M0 and M1 are arguments to this function
729  // the vectors Z0, Z1 : GetDXDAZ()
730  //
731  // The matrix A is calculated from a matrix B as
732  //
733  // A_ij = B_ij / sum_k B_kj
734  //
735  // where k runs over additional indices of B, not present in A.
736  // (underflow and overflow bins, used for efficiency corrections)
737  //
738  // define: Norm_j = sum_k B_kj (data member fSumOverY)
739  //
740  // the derivative of A wrt this input matrix B is given by:
741  //
742  // dA_ij / dB_kj = ( delta_ik - A_ij ) * 1/Norm_j
743  //
744  // The covariance matrix Vxx is:
745  //
746  // Vxx_mn = sum_ijlk [ (dX_m / dA_ij) * (dA_ij / dB_kj) * DB_kj
747  // * (dX_n / dA_lj) * (dA_lj / dB_kj) ]
748  //
749  // where DB_kj is the error on B_kj squared
750  // Simplify the sum over k:
751  //
752  // sum_k [ (dA_ij / dB_kj) * DB_kj * (dA_lj / dB_kj) ]
753  // = sum_k [ ( delta_ik - A_ij ) * 1/Norm_j * DB_kj *
754  // * ( delta_lk - A_lj ) * 1/Norm_j ]
755  // = sum_k [ ( delta_ik*delta_lk - delta_ik*A_lj - delta_lk*A_ij
756  // + A_ij * A_lj ) * DB_kj / Norm_j^2 ]
757  //
758  // introduce normalized errors: Rsq_kj = DB_kj / Norm_j^2
759  // after summing over k:
760  // delta_ik*delta_lk*Rsq_kj -> delta_il*Rsq_ij
761  // delta_ik*A_lj*Rsq_kj -> A_lj*Rsq_ij
762  // delta_lk*A_ij*Rsq_kj -> A_ij*Rsq_lj
763  // A_ij*A_lj*Rsq_kj -> A_ij*A_lj*sum_k(Rsq_kj)
764  //
765  // introduce sum of normalized errors squared: SRsq_j = sum_k(Rsq_kj)
766  //
767  // Note: Rsq_ij is stored as fDAinRelSq (excludes extra indices of B)
768  // and SRsq_j is stored as fDAinColRelSq (sum includes all indices of B)
769  //
770  // Vxx_nm = sum_ijl [ (dX_m / dA_ij) * (dX_n / dA_lj)
771  // (delta_il*Rsq_ij - A_lj*Rsq_ij - A_ij*Rsq_lj + A_ij*A_lj *SRsq_j) ]
772  //
773  // Vxx_nm = sum_j [ F_mj * F_nj * SRsq_j
774  // - sum_j [ G_mj * F_nj ]
775  // - sum_j [ F_mj * G_nj ]
776  // + sum_ij [ (dX_m / dA_ij) * (dX_n / dA_lj) * Rsq_ij ]
777  //
778  // where:
779  // F_mj = sum_i [ (dX_m / dA_ij) * A_ij ]
780  // G_mj = sum_i [ (dX_m / dA_ij) * Rsq_ij ]
781  //
782  // In order to avoid explicitly calculating the 3-dimensional tensor
783  // (dX_m/dA_ij) the sums are evaluated further, using
784  // dX_k / dA_ij = M0_kj * Z0_i - M1_ki * Z1_j
785  //
786  // F_mj = M0_mj * (A# Z0)_j - (M1 A)_mj Z1_j
787  // G_mj = M0_mj * (Rsq# Z0)_j - (M1 Rsq)_mj Z1_j
788  //
789  // and
790  //
791  // sum_ij [ (dX_m/dA_ij) * (dX_n/dA_ij) * Rsq_ij ] =
792  // sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
793  // + sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
794  // - sum_i [ M1_mi * H_ni + M1_ni * H_mi]
795  // where:
796  // H_mi = Z0_i * sum_j [ M0_mj * Z1_j * Rsq_ij ]
797  //
798  // collect all contributions:
799  // Vxx_nm = r0 -r1 -r2 +r3 +r4 -r5 -r6
800  // r0 = sum_j [ F_mj * F_nj * SRsq_j ]
801  // r1 = sum_j [ G_mj * F_nj ]
802  // r2 = sum_j [ F_mj * G_nj ]
803  // r3 = sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
804  // r4 = sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
805  // r5 = sum_i [ M1_mi * H_ni ]
806  // r6 = sum_i [ M1_ni * H_mi ]
807 
808  //======================================================
809  // calculate contributions containing matrices F and G
810  // r0,r1,r2
811  TMatrixDSparse *r=0;
812  if(fDAinColRelSq && fDAinRelSq) {
813  // calculate matrices (M1*A)_mj * Z1_j and (M1*Rsq)_mj * Z1_j
814  TMatrixDSparse *M1A_Z1=MultiplyMSparseMSparse(m_1,fA);
815  ScaleColumnsByVector(M1A_Z1,GetDXDAZ(1));
816  TMatrixDSparse *M1Rsq_Z1=MultiplyMSparseMSparse(m_1,fDAinRelSq);
817  ScaleColumnsByVector(M1Rsq_Z1,GetDXDAZ(1));
818  // calculate vectors A#*Z0 and Rsq#*Z0
819  TMatrixDSparse *AtZ0 = MultiplyMSparseTranspMSparse(fA,GetDXDAZ(0));
820  TMatrixDSparse *RsqZ0=
821  MultiplyMSparseTranspMSparse(fDAinRelSq,GetDXDAZ(0));
822  //calculate matrix F
823  // F_mj = M0_mj * (A# Z0)_j - (M1 A)_mj Z1_j
824  TMatrixDSparse *F=new TMatrixDSparse(*m_0);
825  ScaleColumnsByVector(F,AtZ0);
826  AddMSparse(F,-1.0,M1A_Z1);
827  //calculate matrix G
828  // G_mj = M0_mj * (Rsq# Z0)_j - (M1 Rsq)_mj Z1_j
829  TMatrixDSparse *G=new TMatrixDSparse(*m_0);
830  ScaleColumnsByVector(G,RsqZ0);
831  AddMSparse(G,-1.0,M1Rsq_Z1);
832  DeleteMatrix(&M1A_Z1);
833  DeleteMatrix(&M1Rsq_Z1);
834  DeleteMatrix(&AtZ0);
835  DeleteMatrix(&RsqZ0);
836  // r0 = sum_j [ F_mj * F_nj * SRsq_j ]
837  r=MultiplyMSparseMSparseTranspVector(F,F,fDAinColRelSq);
838  // r1 = sum_j [ G_mj * F_nj ]
839  TMatrixDSparse *r1=MultiplyMSparseMSparseTranspVector(F,G,0);
840  // r2 = sum_j [ F_mj * G_nj ]
841  TMatrixDSparse *r2=MultiplyMSparseMSparseTranspVector(G,F,0);
842  // r = r0-r1-r2
843  AddMSparse(r,-1.0,r1);
844  AddMSparse(r,-1.0,r2);
845  DeleteMatrix(&r1);
846  DeleteMatrix(&r2);
847  DeleteMatrix(&F);
848  DeleteMatrix(&G);
849  }
850  //======================================================
851  // calculate contribution
852  // sum_ij [ (dX_m/dA_ij) * (dX_n/dA_ij) * Rsq_ij ]
853  // (r3,r4,r5,r6)
854  if(fDAinRelSq) {
855  // (Z0_i)^2
856  TMatrixDSparse Z0sq(*GetDXDAZ(0));
857  const Int_t *Z0sq_rows=Z0sq.GetRowIndexArray();
858  Double_t *Z0sq_data=Z0sq.GetMatrixArray();
859  for(int index=0;index<Z0sq_rows[Z0sq.GetNrows()];index++) {
860  Z0sq_data[index] *= Z0sq_data[index];
861  }
862  // Z0sqRsq = sum_i (Z_i)^2 * Rsq_ij
863  TMatrixDSparse *Z0sqRsq=MultiplyMSparseTranspMSparse(fDAinRelSq,&Z0sq);
864  // r3 = sum_j [ M0_mj * M0_nj * [ sum_i (Z0_i)^2 * Rsq_ij ] ]
865  TMatrixDSparse *r3=MultiplyMSparseMSparseTranspVector(m_0,m_0,Z0sqRsq);
866  DeleteMatrix(&Z0sqRsq);
867 
868  // (Z1_j)^2
869  TMatrixDSparse Z1sq(*GetDXDAZ(1));
870  const Int_t *Z1sq_rows=Z1sq.GetRowIndexArray();
871  Double_t *Z1sq_data=Z1sq.GetMatrixArray();
872  for(int index=0;index<Z1sq_rows[Z1sq.GetNrows()];index++) {
873  Z1sq_data[index] *= Z1sq_data[index];
874  }
875  // Z1sqRsq = sum_j (Z1_j)^2 * Rsq_ij ]
876  TMatrixDSparse *Z1sqRsq=MultiplyMSparseMSparse(fDAinRelSq,&Z1sq);
877  // r4 = sum_i [ M1_mi * M1_ni * [ sum_j (Z1_j)^2 * Rsq_ij ] ]
878  TMatrixDSparse *r4=MultiplyMSparseMSparseTranspVector(m_1,m_1,Z1sqRsq);
879  DeleteMatrix(&Z1sqRsq);
880 
881  // sum_j [ M0_mj * Z1_j * Rsq_ij ]
882  TMatrixDSparse *H=MultiplyMSparseMSparseTranspVector
883  (m_0,fDAinRelSq,GetDXDAZ(1));
884  // H_mi = Z0_i * sum_j [ M0_mj * Z1_j * Rsq_ij ]
885  ScaleColumnsByVector(H,GetDXDAZ(0));
886  // r5 = sum_i [ M1_mi * H_ni ]
887  TMatrixDSparse *r5=MultiplyMSparseMSparseTranspVector(m_1,H,0);
888  // r6 = sum_i [ H_mi * M1_ni ]
889  TMatrixDSparse *r6=MultiplyMSparseMSparseTranspVector(H,m_1,0);
890  DeleteMatrix(&H);
891  // r = r0 -r1 -r2 +r3 +r4 +r5 +r6
892  if(r) {
893  AddMSparse(r,1.0,r3);
894  DeleteMatrix(&r3);
895  } else {
896  r=r3;
897  r3=0;
898  }
899  AddMSparse(r,1.0,r4);
900  AddMSparse(r,-1.0,r5);
901  AddMSparse(r,-1.0,r6);
902  DeleteMatrix(&r4);
903  DeleteMatrix(&r5);
904  DeleteMatrix(&r6);
905  }
906  return r;
907 }
908 
910 (const TMatrixDSparse *m1,const TMatrixDSparse *m2,const TMatrixDSparse *dsys)
911 {
912  // propagate correlated systematic shift to output vector
913  // m1,m2 : coefficients for propagating the errors
914  // dsys : matrix of correlated shifts from this source
915 
916  // delta_m =
917  // sum{i,j} {
918  // ((*m1)(m,j) * (*fVYAx)(i) - (*m2)(m,i) * (*fX)(j))*dsys(i,j) }
919  // = sum_j (*m1)(m,j) sum_i dsys(i,j) * (*fVYAx)(i)
920  // - sum_i (*m2)(m,i) sum_j dsys(i,j) * (*fX)(j)
921 
922  TMatrixDSparse *dsysT_VYAx = MultiplyMSparseTranspMSparse(dsys,GetDXDAZ(0));
923  TMatrixDSparse *delta = MultiplyMSparseMSparse(m1,dsysT_VYAx);
924  DeleteMatrix(&dsysT_VYAx);
925  TMatrixDSparse *dsys_X = MultiplyMSparseMSparse(dsys,GetDXDAZ(1));
926  TMatrixDSparse *delta2 = MultiplyMSparseMSparse(m2,dsys_X);
927  DeleteMatrix(&dsys_X);
928  AddMSparse(delta,-1.0,delta2);
929  DeleteMatrix(&delta2);
930  return delta;
931 }
932 
934 {
935  // set uncertainty on tau
936  fDtau=delta_tau;
938 }
939 
940 Bool_t TUnfoldSys::GetDeltaSysSource(TH1 *hist_delta,const char *name,
941  const Int_t *binMap)
942 {
943  // calculate systematic shift from a given source
944  // ematrix: output
945  // source: name of the error source
946  // binMap: see method GetEmatrix()
947  PrepareSysError();
948  const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
949  const TMatrixDSparse *delta=0;
950  if(named_emat) {
951  delta=(TMatrixDSparse *)named_emat->Value();
952  }
953  VectorMapToHist(hist_delta,delta,binMap);
954  return delta !=0;
955 }
956 
958 (TH1 *hist_delta,const char *source,const Int_t *binMap)
959 {
960  // get correlated shift induced by a background source
961  // delta: output shift vector histogram
962  // source: name of background source
963  // binMap: see method GetEmatrix()
964  // see PrepareSysError()
965  PrepareSysError();
966  const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(source);
967  TMatrixDSparse *dx=0;
968  if(named_err) {
969  const TMatrixD *dy=(TMatrixD *)named_err->Value();
970  dx=MultiplyMSparseM(GetDXDY(),dy);
971  }
972  VectorMapToHist(hist_delta,dx,binMap);
973  if(dx!=0) {
974  DeleteMatrix(&dx);
975  return kTRUE;
976  }
977  return kFALSE;
978 }
979 
980 Bool_t TUnfoldSys::GetDeltaSysTau(TH1 *hist_delta,const Int_t *binMap)
981 {
982  // calculate systematic shift from tau variation
983  // ematrix: output
984  // binMap: see method GetEmatrix()
985  PrepareSysError();
986  VectorMapToHist(hist_delta,fDeltaSysTau,binMap);
987  return fDeltaSysTau !=0;
988 }
989 
991 (TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
992 {
993  // calculate systematic shift from a given source
994  // ematrix: output
995  // source: name of the error source
996  // binMap: see method GetEmatrix()
997  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
998  PrepareSysError();
999  const TPair *named_emat=(const TPair *)fDeltaCorrX->FindObject(name);
1000  TMatrixDSparse *emat=0;
1001  if(named_emat) {
1002  const TMatrixDSparse *delta=(TMatrixDSparse *)named_emat->Value();
1003  emat=MultiplyMSparseMSparseTranspVector(delta,delta,0);
1004  }
1005  ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1006  DeleteMatrix(&emat);
1007 }
1008 
1010 (TH2 *ematrix,const char *name,const Int_t *binMap,Bool_t clearEmat)
1011 {
1012  // calculate systematic shift from a given background scale error
1013  // ematrix: output
1014  // source: name of the error source
1015  // binMap: see method GetEmatrix()
1016  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1017  PrepareSysError();
1018  const TPair *named_err=(const TPair *)fBgrErrScaleIn->FindObject(name);
1019  TMatrixDSparse *emat=0;
1020  if(named_err) {
1021  const TMatrixD *dy=(TMatrixD *)named_err->Value();
1022  TMatrixDSparse *dx=MultiplyMSparseM(GetDXDY(),dy);
1023  emat=MultiplyMSparseMSparseTranspVector(dx,dx,0);
1024  DeleteMatrix(&dx);
1025  }
1026  ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1027  DeleteMatrix(&emat);
1028 }
1029 
1031 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1032 {
1033  // calculate error matrix from error in regularisation parameter
1034  // ematrix: output
1035  // binMap: see method GetEmatrix()
1036  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1037  PrepareSysError();
1038  TMatrixDSparse *emat=0;
1039  if(fDeltaSysTau) {
1040  emat=MultiplyMSparseMSparseTranspVector(fDeltaSysTau,fDeltaSysTau,0);
1041  }
1042  ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1043  DeleteMatrix(&emat);
1044 }
1045 
1047 (TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1048 {
1049  // calculate error matrix from error in input vector alone
1050  // ematrix: output
1051  // binMap: see method GetEmatrix()
1052  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1053  GetEmatrixFromVyy(fVyyData,ematrix,binMap,clearEmat);
1054 }
1055 
1057 (TH2 *ematrix,const char *source,const Int_t *binMap,Bool_t clearEmat)
1058 {
1059  // calculate error matrix contribution originating from uncorrelated errors
1060  // of one background source
1061  // ematrix: output
1062  // source: name of the error source
1063  // binMap: see method GetEmatrix()
1064  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1065  const TPair *named_err=(const TPair *)fBgrErrUncorrInSq->FindObject(source);
1066  TMatrixDSparse *emat=0;
1067  if(named_err) {
1068  TMatrixD const *dySquared=(TMatrixD const *)named_err->Value();
1069  emat=MultiplyMSparseMSparseTranspVector(GetDXDY(),GetDXDY(),dySquared);
1070  }
1071  ErrorMatrixToHist(ematrix,emat,binMap,clearEmat);
1072  DeleteMatrix(&emat);
1073 }
1074 
1076 (const TMatrixDSparse *vyy,TH2 *ematrix,const Int_t *binMap,Bool_t clearEmat)
1077 {
1078  // propagate error matrix vyy to the result
1079  // vyy: error matrix on input data fY
1080  // ematrix: output
1081  // binMap: see method GetEmatrix()
1082  // clearEmat: set kTRUE to clear the histogram prior to adding the errors
1083  PrepareSysError();
1084  TMatrixDSparse *em=0;
1085  if(vyy) {
1086  TMatrixDSparse *dxdyVyy=MultiplyMSparseMSparse(GetDXDY(),vyy);
1087  em=MultiplyMSparseMSparseTranspVector(dxdyVyy,GetDXDY(),0);
1088  DeleteMatrix(&dxdyVyy);
1089  }
1090  ErrorMatrixToHist(ematrix,em,binMap,clearEmat);
1091  DeleteMatrix(&em);
1092 }
1093 
1094 void TUnfoldSys::GetEmatrixTotal(TH2 *ematrix,const Int_t *binMap)
1095 {
1096  // get total error including statistical error
1097  // ematrix: output
1098  // binMap: see method GetEmatrix()
1099  GetEmatrix(ematrix,binMap); // (stat)+(d)+(e)
1100  GetEmatrixSysUncorr(ematrix,binMap,kFALSE); // (a)
1101  TMapIter sysErrPtr(fDeltaCorrX);
1102  const TObject *key;
1103 
1104  for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1105  GetEmatrixSysSource(ematrix,
1106  ((const TObjString *)key)->GetString(),
1107  binMap,kFALSE); // (b)
1108  }
1109  GetEmatrixSysTau(ematrix,binMap,kFALSE); // (c)
1110 }
1111 
1113 {
1114  PrepareSysError();
1115 
1116  // errors from input vector and from background subtraction
1117  TMatrixDSparse *emat_sum=new TMatrixDSparse(*fVyy);
1118 
1119  // uncorrelated systematic error
1120  AddMSparse(emat_sum,1.0,fEmatUncorrAx);
1121  TMapIter sysErrPtr(fDeltaCorrAx);
1122  const TObject *key;
1123 
1124  // correlated systematic errors
1125  for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1126  const TMatrixDSparse *delta=(TMatrixDSparse *)
1127 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
1128  ((const TPair *)*sysErrPtr)->Value()
1129 #else
1130  fDeltaCorrAx->GetValue(((const TObjString *)key)
1131  ->GetString())
1132 #endif
1133  ;
1135  AddMSparse(emat_sum,1.0,emat);
1136  DeleteMatrix(&emat);
1137  }
1138  // error on tau
1139  if(fDeltaSysTau) {
1141  TMatrixDSparse *emat_tau=
1142  MultiplyMSparseMSparseTranspVector(Adx_tau,Adx_tau,0);
1143  DeleteMatrix(&Adx_tau);
1144  AddMSparse(emat_sum,1.0,emat_tau);
1145  DeleteMatrix(&emat_tau);
1146  }
1147  return emat_sum;
1148 }
1149 
1151 {
1152  PrepareSysError();
1153 
1154  // errors from input vector and from background subtraction
1155  TMatrixDSparse *emat_sum=new TMatrixDSparse(*GetVxx());
1156 
1157  // uncorrelated systematic error
1158  AddMSparse(emat_sum,1.0,fEmatUncorrX);
1159  TMapIter sysErrPtr(fDeltaCorrX);
1160  const TObject *key;
1161 
1162  // correlated systematic errors
1163  for(key=sysErrPtr.Next();key;key=sysErrPtr.Next()) {
1164  const TMatrixDSparse *delta=(TMatrixDSparse *)
1165 #if ROOT_VERSION_CODE >= ROOT_VERSION(5,20,00)
1166  ((const TPair *)*sysErrPtr)->Value()
1167 #else
1168  fDeltaCorrX->GetValue(((const TObjString *)key)
1169  ->GetString())
1170 #endif
1171  ;
1173  AddMSparse(emat_sum,1.0,emat);
1174  DeleteMatrix(&emat);
1175  }
1176  // error on tau
1177  if(fDeltaSysTau) {
1178  TMatrixDSparse *emat_tau=
1180  AddMSparse(emat_sum,1.0,emat_tau);
1181  DeleteMatrix(&emat_tau);
1182  }
1183  return emat_sum;
1184 }
1185 
1186 
1188 {
1189  // calculate total chi**2 including systematic errors
1190 
1192 
1193  Int_t rank=0;
1194  TMatrixDSparse *v=InvertMSparseSymmPos(emat_sum,&rank);
1195  TMatrixD dy(*fY, TMatrixD::kMinus, *GetAx());
1196  TMatrixDSparse *vdy=MultiplyMSparseM(v,&dy);
1197  DeleteMatrix(&v);
1198  Double_t r=0.0;
1199  const Int_t *vdy_rows=vdy->GetRowIndexArray();
1200  const Double_t *vdy_data=vdy->GetMatrixArray();
1201  for(Int_t i=0;i<vdy->GetNrows();i++) {
1202  if(vdy_rows[i+1]>vdy_rows[i]) {
1203  r += vdy_data[vdy_rows[i]] * dy(i,0);
1204  }
1205  }
1206  DeleteMatrix(&vdy);
1207  DeleteMatrix(&emat_sum);
1208  return r;
1209 }
1210 
1211 void TUnfoldSys::GetRhoItotal(TH1 *rhoi,const Int_t *binMap,TH2 *invEmat)
1212 {
1213  // get global correlation coefficients including systematic,statistical,background,tau errors
1214  // rhoi: output histogram
1215  // binMap: for each global bin, indicate in which histogram bin
1216  // to store its content
1217  // invEmat: output histogram for inverse of error matrix
1218  // (pointer may zero if inverse is not requested)
1219  ClearHistogram(rhoi,-1.);
1221  GetRhoIFromMatrix(rhoi,emat_sum,binMap,invEmat);
1222 
1223  DeleteMatrix(&emat_sum);
1224 }
1225 
1228 {
1229  // scale columns of m by the corresponding rows of v
1230  // input:
1231  // m: pointer to sparse matrix of dimension NxM
1232  // v: pointer to matrix of dimension Mx1
1233  if((m->GetNcols() != v->GetNrows())||(v->GetNcols()!=1)) {
1234  Fatal("ScaleColumnsByVector error",
1235  "matrix cols/vector rows %d!=%d OR vector cols %d !=1\n",
1236  m->GetNcols(),v->GetNrows(),v->GetNcols());
1237  }
1238  const Int_t *rows_m=m->GetRowIndexArray();
1239  const Int_t *cols_m=m->GetColIndexArray();
1240  Double_t *data_m=m->GetMatrixArray();
1241  const TMatrixDSparse *v_sparse=dynamic_cast<const TMatrixDSparse *>(v);
1242  if(v_sparse) {
1243  const Int_t *rows_v=v_sparse->GetRowIndexArray();
1244  const Double_t *data_v=v_sparse->GetMatrixArray();
1245  for(Int_t i=0;i<m->GetNrows();i++) {
1246  for(Int_t index_m=rows_m[i];index_m<rows_m[i+1];index_m++) {
1247  Int_t j=cols_m[index_m];
1248  Int_t index_v=rows_v[j];
1249  if(index_v<rows_v[j+1]) {
1250  data_m[index_m] *= data_v[index_v];
1251  } else {
1252  data_m[index_m] =0.0;
1253  }
1254  }
1255  }
1256  } else {
1257  for(Int_t i=0;i<m->GetNrows();i++) {
1258  for(Int_t index=rows_m[i];index<rows_m[i+1];index++) {
1259  data_m[index] *= (*v)(cols_m[index],0);
1260  }
1261  }
1262  }
1263 }
1264 
1266 (TH1 *hist_delta,const TMatrixDSparse *delta,const Int_t *binMap)
1267 {
1268  // sum over bins of *delta, as defined in binMap,fXToHist
1269  // hist_delta: histogram to return summed vector
1270  // delta: vector to sum and remap
1271  Int_t nbin=hist_delta->GetNbinsX();
1272  Double_t *c=new Double_t[nbin+2];
1273  for(Int_t i=0;i<nbin+2;i++) {
1274  c[i]=0.0;
1275  }
1276  if(delta) {
1277  Int_t binMapSize = fHistToX.GetSize();
1278  const Double_t *delta_data=delta->GetMatrixArray();
1279  const Int_t *delta_rows=delta->GetRowIndexArray();
1280  for(Int_t i=0;i<binMapSize;i++) {
1281  Int_t destBinI=binMap ? binMap[i] : i;
1282  Int_t srcBinI=fHistToX[i];
1283  if((destBinI>=0)&&(destBinI<nbin+2)&&(srcBinI>=0)) {
1284  Int_t index=delta_rows[srcBinI];
1285  if(index<delta_rows[srcBinI+1]) {
1286  c[destBinI]+=delta_data[index];
1287  }
1288  }
1289  }
1290  }
1291  for(Int_t i=0;i<nbin+2;i++) {
1292  hist_delta->SetBinContent(i,c[i]);
1293  hist_delta->SetBinError(i,0.0);
1294  }
1295  delete[] c;
1296 }
void InitTUnfoldSys(void)
Definition: TUnfoldSys.cxx:551
virtual Int_t GetEntries() const
Definition: TCollection.h:92
TMatrixDSparse * fEmatUncorrX
Definition: TUnfoldSys.h:65
static void DeleteMatrix(TMatrixD **m)
Definition: TUnfold.cxx:333
TMatrixDSparse * CreateSparseMatrix(Int_t nrow, Int_t ncol, Int_t nele, Int_t *row, Int_t *col, Double_t *data) const
Definition: TUnfold.cxx:698
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
virtual void ClearResults(void)
Definition: TUnfoldSys.cxx:608
TMatrixDSparse * InvertMSparseSymmPos(const TMatrixDSparse *A, Int_t *rank) const
Definition: TUnfold.cxx:1064
Double_t GetChi2Sys(void)
void Fatal(const char *location, const char *msgfmt,...)
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Collectable string class.
Definition: TObjString.h:32
TMap * fBgrErrScaleIn
Definition: TUnfoldSys.h:61
virtual const Element * GetMatrixArray() const
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:940
TMap * fBgrIn
Definition: TUnfoldSys.h:59
void VectorMapToHist(TH1 *hist_delta, const TMatrixDSparse *delta, const Int_t *binMap)
virtual TMatrixDSparse * PrepareCorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixDSparse *dsys)
Definition: TUnfoldSys.cxx:910
virtual void SetOwner(Bool_t enable=kTRUE)
Set whether this collection is the owner (enable==true) of its content.
TMatrixDSparse * fDAinRelSq
Definition: TUnfoldSys.h:55
virtual void PrepareSysError(void)
Definition: TUnfoldSys.cxx:619
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:52
#define H(x, y, z)
TUnfold is used to decompose a measurement y into several sources x given the measurement uncertainti...
Definition: TUnfoldSys.h:51
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
const Bool_t kFALSE
Definition: Rtypes.h:92
void SetTauError(Double_t delta_tau)
Definition: TUnfoldSys.cxx:933
Double_t fDtau
Definition: TUnfoldSys.h:62
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)
Definition: TUnfoldSys.cxx:410
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
#define ROOT_VERSION_CODE
Definition: RVersion.h:21
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
TMap * fSysIn
Definition: TUnfoldSys.h:58
#define G(x, y, z)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
const TMatrixDSparse * GetDXDtauSquared(void) const
Definition: TUnfold.h:173
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:946
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse *a, const TMatrixD *b) const
Definition: TUnfold.cxx:858
void GetEmatrixSysUncorr(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Definition: TUnfoldSys.cxx:704
void ClearHistogram(TH1 *h, Double_t x=0.) const
Definition: TUnfold.cxx:3394
unsigned int r3[N_CITIES]
Definition: simanTSP.cxx:323
const int ny
Definition: kalman.C:17
EHistMap
Definition: TUnfold.h:188
EConstraint
Definition: TUnfold.h:103
TMatrixDSparse * fDeltaSysTau
Definition: TUnfoldSys.h:69
void GetEmatrixSysBackgroundScale(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
void Info(const char *location, const char *msgfmt,...)
TObject * GetValue(const char *keyname) const
Returns a pointer to the value associated with keyname as name of the key.
Definition: TMap.cxx:234
void DoBackgroundSubtraction(void)
Definition: TUnfoldSys.cxx:300
const TMatrixDSparse * GetVxx(void) const
Definition: TUnfold.h:177
TMatrixDSparse * GetSummedErrorMatrixXX(void)
virtual void SetBinError(Int_t bin, Double_t error)
see convention for numbering bins in TH1::GetBin
Definition: TH1.cxx:8528
TMatrixDSparse * GetSummedErrorMatrixYY(void)
void Error(const char *location, const char *msgfmt,...)
TMatrixTSparse< Double_t > TMatrixDSparse
ERegMode
Definition: TUnfold.h:107
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:90
TObject * FindObject(const char *keyname) const
Check if a (key,value) pair exists with keyname as name of the key.
Definition: TMap.cxx:213
void GetEmatrix(TH2 *ematrix, const Int_t *binMap=0) const
Definition: TUnfold.cxx:3175
TMatrixT< Double_t > TMatrixD
Definition: TMatrixDfwd.h:24
Iterator of map.
Definition: TMap.h:148
TMap * fBgrErrUncorrInSq
Definition: TUnfoldSys.h:60
#define F(x, y, z)
virtual void SetOwnerKeyValue(Bool_t ownkeys=kTRUE, Bool_t ownvals=kTRUE)
Set ownership for keys and values.
Definition: TMap.cxx:350
void GetEmatrixSysSource(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Definition: TUnfoldSys.cxx:991
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)
Definition: TUnfold.cxx:2208
TMap * fDeltaCorrX
Definition: TUnfoldSys.h:67
const TMatrixDSparse * GetAx(void) const
Definition: TUnfold.h:174
Int_t GetNy(void) const
Definition: TUnfold.h:165
TString GetString() const
Definition: TObjString.h:50
TMap * fDeltaCorrAx
Definition: TUnfoldSys.h:68
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual const Int_t * GetRowIndexArray() const
Service class for 2-Dim histogram classes.
Definition: TH2.h:36
SVector< double, 2 > v
Definition: Dict.h:5
virtual TMatrixDSparse * PrepareUncorrEmat(const TMatrixDSparse *m1, const TMatrixDSparse *m2)
Definition: TUnfoldSys.cxx:717
void GetEmatrixTotal(TH2 *ematrix, const Int_t *binMap=0)
unsigned int r1[N_CITIES]
Definition: simanTSP.cxx:321
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:8543
TObject * Value() const
Definition: TMap.h:125
const TMatrixDSparse * GetVyyInv(void) const
Definition: TUnfold.h:179
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:980
virtual ~TUnfoldSys(void)
Definition: TUnfoldSys.cxx:592
TMatrixD * fY
Definition: TUnfold.h:118
Double_t fTauSquared
Definition: TUnfold.h:120
TMatrixDSparse * fVyy
Definition: TUnfold.h:117
void SubtractBackground(const TH1 *hist_bgr, const char *name, Double_t scale=1.0, Double_t scale_error=0.0)
Definition: TUnfoldSys.cxx:441
TMatrixDSparse * fVyyData
Definition: TUnfoldSys.h:64
Class used by TMap to store (key,value) pairs.
Definition: TMap.h:106
void GetEmatrixSysTau(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
double Double_t
Definition: RtypesCore.h:55
virtual void ClearResults(void)
Definition: TUnfold.cxx:345
TMatrixD * fDAinColRelSq
Definition: TUnfoldSys.h:56
TMap implements an associative array of (key,value) pairs using a THashTable for efficient retrieval ...
Definition: TMap.h:44
Int_t GetNcols() const
Definition: TMatrixTBase.h:137
TUnfoldSys(void)
The TH1 histogram class.
Definition: TH1.h:80
TMatrixDSparse * fEmatUncorrAx
Definition: TUnfoldSys.h:66
void AddSysError(const TH2 *sysError, const char *name, EHistMap histmap, ESysErrMode mode)
Definition: TUnfoldSys.cxx:217
#define ROOT_VERSION(a, b, c)
Definition: RVersion.h:20
Mother of all ROOT objects.
Definition: TObject.h:58
virtual const Int_t * GetColIndexArray() const
TMatrixDSparse * fA
Definition: TUnfold.h:115
void AddMSparse(TMatrixDSparse *dest, Double_t f, const TMatrixDSparse *src) const
Definition: TUnfold.cxx:1001
ClassImp(TUnfoldSys) TUnfoldSys
Definition: TUnfoldSys.cxx:133
TUnfold is used to decompose a measurement y into several sources x given the measurement uncertainti...
Definition: TUnfold.h:99
TMatrixDSparse * MultiplyMSparseMSparse(const TMatrixDSparse *a, const TMatrixDSparse *b) const
Definition: TUnfold.cxx:711
TObject * Next()
Returns the next key from a map.
Definition: TMap.cxx:531
void ScaleColumnsByVector(TMatrixDSparse *m, const TMatrixTBase< Double_t > *v) const
TMatrixD * fYData
Definition: TUnfoldSys.h:63
const TMatrixDSparse * GetDXDAM(int i) const
Definition: TUnfold.h:171
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
void GetRhoItotal(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0)
virtual Double_t GetBinError(Int_t bin) const
Return value of error associated to bin number bin.
Definition: TH1.cxx:8395
const Bool_t kTRUE
Definition: Rtypes.h:91
void GetEmatrixFromVyy(const TMatrixDSparse *vyy, TH2 *ematrix, const Int_t *binMap, Bool_t clearEmat)
Double_t GetRhoIFromMatrix(TH1 *rhoi, const TMatrixDSparse *eOrig, const Int_t *binMap, TH2 *invEmat) const
Definition: TUnfold.cxx:3275
virtual TObject * FindObject(const char *name) const
Must be redefined in derived classes.
Definition: TObject.cxx:379
const Int_t n
Definition: legend1.C:16
Bool_t GetDeltaSysBackgroundScale(TH1 *delta, const char *source, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:958
unsigned int r2[N_CITIES]
Definition: simanTSP.cxx:322
int CompareTo(const char *cs, ECaseCompare cmp=kExact) const
Compare a string to char *cs2.
Definition: TString.cxx:385
void Clear(Option_t *option="")
Remove all (key,value) pairs from the map.
Definition: TMap.cxx:95
TMatrixD * fAoutside
Definition: TUnfoldSys.h:57
const char * Value
Definition: TXMLSetup.cxx:73
TMatrixDSparse * MultiplyMSparseMSparseTranspVector(const TMatrixDSparse *m1, const TMatrixDSparse *m2, const TMatrixTBase< Double_t > *v) const
Definition: TUnfold.cxx:911
void GetBackground(TH1 *bgr, const char *bgrSource=0, const Int_t *binMap=0, Int_t includeError=3, Bool_t clearHist=kTRUE) const
Definition: TUnfoldSys.cxx:480
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:904