Logo ROOT   6.14/05
Reference Guide
TUnfoldDensity.cxx
Go to the documentation of this file.
1 // @(#)root/unfold:$Id$
2 // Authors: Stefan Schmitt, Amnon Harel DESY and CERN, 11/08/11
3 
4 /** \class TUnfoldDensity
5 \ingroup Unfold
6 An algorithm to unfold distributions from detector to truth level
7 
8 TUnfoldDensity is used to decompose a measurement y into several sources x,
9 given the measurement uncertainties, background b and a matrix of migrations A.
10 The method can be applied to a large number of problems,
11 where the measured distribution y is a linear superposition
12 of several Monte Carlo shapes. Beyond such a simple template fit,
13 TUnfoldDensity has an adjustable regularisation term and also supports an
14 optional constraint on the total number of events.
15 Background sources can be specified, with a normalisation constant and
16 normalisation uncertainty. In addition, variants of the response
17 matrix may be specified, these are taken to determine systematic
18 uncertainties. Complex, multidimensional arrangements of signal and
19 background bins are managed with the help of the class TUnfoldBinning.
20 
21 If you use this software, please consider the following citation
22 
23 <b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
24 
25 Detailed documentation and updates are available on
26 http://www.desy.de/~sschmitt
27 
28 ### Brief recipe to use TUnfoldSys:
29 
30  - Set up binning schemes for the truth and measured
31 distributions. The binning schemes may be coded in the XML language,
32 for reading use TUnfoldBinningXML.
33  - A matrix (truth,reconstructed) is given as a two-dimensional histogram
34  as argument to the constructor of TUnfold
35  - A vector of measurements is given as one-dimensional histogram using
36  the SetInput() method
37  - Repeated calls to SubtractBackground() to specify background sources
38  - Repeated calls to AddSysError() to specify systematic uncertainties
39  - The unfolding is performed
40 
41  - either once with a fixed parameter tau, method DoUnfold(tau)
42  - or multiple times in a scan to determine the best choice of tau,
43  method ScanLCurve()
44  - or multiple times in a scan to determine the best choice of tau,
45  method ScanTau()
46 
47  - Unfolding results are retrieved using various GetXXX() methods
48 
49 A detailed documentation of the various GetXXX() methods to control
50 systematic uncertainties is given with the method TUnfoldSys.
51 
52 ### Why to use complex binning schemes
53 
54 in literature on unfolding, the "standard" test case is a
55 one-dimensional distribution without underflow or overflow bins.
56 The migration matrix is almost diagonal.
57 
58 <b>This "standard" case is rarely realized for real problems.</b>
59 
60 Often one has to deal with multi-dimensional distributions.
61 In addition, there are underflow and overflow bins
62 or other background bins, possibly determined with the help of auxiliary
63 measurements.
64 
65 In TUnfoldDensity, such complex binning schemes are handled with the help
66 of the class TUnfoldBinning. For both the measurement and the truth
67 there is a tree structure. The tree nodes may correspond to single
68 bins (e.g. nuisance parameters) or may hold multi-dimensional distributions.
69 
70 For example, the "measurement" tree could have two leaves, one for
71 the primary distribution and one for auxiliary measurements.
72 Similarly, the "truth" tree could have two leaves, one for the
73 signal and one for the background.
74 Each of the leaves may then have a multi-dimensional distribution.
75 
76 The class TUnfoldBinning takes care to map all bins of the
77 "measurement" to a one-dimensional vector y.
78 Similarly, the "truth" bins are mapped to the vector x.
79 
80 ### How to choose the regularisation settings
81 
82 In TUnfoldDensity, two methods are implemented to determine tau**2
83 
84  1. ScanLcurve() locate the tau where the L-curve plot has a "kink"
85  this function is implemented in the TUnfold class
86  2. ScanTau() finds the solution such that some variable
87  (e.g. global correlation coefficient) is minimized.
88  This function is implemented in the TUnfoldDensity class
89 
90 Each of the algorithms has its own advantages and disadvantages.
91 The algorithm (1) does not work if the input data are too similar to the
92 MC prediction. Typical no-go cases of the L-curve scan are:
93 
94  - the number of measurements is too small (e.g. ny=nx)
95  - the input data have no statistical fluctuations
96  [identical MC events are used to fill the matrix of migrations
97  and the vector y for a "closure test"]
98 
99 The algorithm (2) only works if the variable does have a real minimum
100 as a function of tau. If global correlations are minimized, the situation
101 is as follows:
102 The matrix of migration typically introduces negative correlations.
103 The area constraint introduces some positive correlation.
104 Regularisation on the "size" introduces no correlation.
105 Regularisation on 1st or 2nd derivatives adds positive correlations.
106 
107 For these reasons, "size" regularisation does not work well with
108 the tau-scan: the higher tau, the smaller rho, but there is no minimum.
109 As a result, large values of tau (too strong regularisation) are found.
110 In contrast, the tau-scan is expected to work better with 1st or 2nd
111 derivative regularisation, because at some point the negative
112 correlations from migrations are approximately cancelled by the
113 positive correlations from the regularisation conditions.
114 
115 whichever algorithm is used, the output has to be checked:
116 
117  1. The L-curve should have approximate L-shape
118  and the final choice of tau should not be at the very edge of the
119  scanned region
120  2. The scan result should have a well-defined minimum and the
121  final choice of tau should sit right in the minimum
122 
123 
124 --------------------------------------------------------------------------------
125  This file is part of TUnfold.
126 
127  TUnfold is free software: you can redistribute it and/or modify
128  it under the terms of the GNU General Public License as published by
129  the Free Software Foundation, either version 3 of the License, or
130  (at your option) any later version.
131 
132  TUnfold is distributed in the hope that it will be useful,
133  but WITHOUT ANY WARRANTY; without even the implied warranty of
134  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
135  GNU General Public License for more details.
136 
137  You should have received a copy of the GNU General Public License
138  along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
139 
140 <b>Version 17.6, with updated doxygen comments and bug-fixes in TUnfoldBinning</b>
141 
142 #### History:
143  - Version 17.5, bug fix in TUnfold also corrects GetEmatrixSysUncorr()
144  - Version 17.4, in parallel to changes in TUnfoldBinning
145  - Version 17.3, in parallel to changes in TUnfoldBinning
146  - Version 17.2, with new options 'N' and 'c' for axis regularisation steering
147  - Version 17.1, add scan type RhoSquare, small bug fixes with useAxisBinning
148  - Version 17.0, support for density regularisation, complex binning schemes, tau scan
149 */
150 
151 #include "TUnfoldDensity.h"
152 #include <TMath.h>
153 #include <TVectorD.h>
154 #include <TObjString.h>
155 #include <iostream>
156 #include <map>
157 
158 //#define DEBUG
159 
160 #ifdef DEBUG
161 using namespace std;
162 #endif
163 
165 
166 TUnfoldDensity::~TUnfoldDensity(void)
167 {
168  // clean up
169  if(fOwnedOutputBins) delete fOwnedOutputBins;
170  if(fOwnedInputBins) delete fOwnedInputBins;
171  if(fRegularisationConditions) delete fRegularisationConditions;
172 }
173 
174 ////////////////////////////////////////////////////////////////////////////////
175 /// Only for use by root streamer or derived classes.
176 
178 {
179  fConstOutputBins=0;
180  fConstInputBins=0;
181  fOwnedOutputBins=0;
182  fOwnedInputBins=0;
183  fRegularisationConditions=0;
184 }
185 
186 ////////////////////////////////////////////////////////////////////////////////
187 /// Eet up response matrix A, uncorrelated uncertainties of A,
188 /// regularisation scheme and binning schemes.
189 ///
190 /// \param[in] hist_A matrix that describes the migrations
191 /// \param[in] histmap mapping of the histogram axes to the unfolding output
192 /// \param[in] regmode (default=kRegModeSize) global regularisation mode
193 /// \param[in] constraint (default=kEConstraintArea) type of constraint
194 /// \param[in] densityMode (default=kDensityModeBinWidthAndUser)
195 /// regularisation scale factors to construct the matrix L
196 /// \param[in] outputBins (default=0) binning scheme for truth (unfolding output)
197 /// \param[in] inputBins (default=0) binning scheme for measurement (unfolding
198 /// input)
199 /// \param[in] regularisationDistribution (default=0) selection of
200 /// regularized distribution
201 /// \param[in] regularisationAxisSteering (default=0) detailed
202 /// regularisation steering for selected distribution
203 ///
204 /// The parameters <b>hist_A, histmap, constraint</b> are
205 /// explained with the TUnfoldSys constructor.
206 ///
207 /// The parameters <b>outputBins,inputBins</b> set the binning
208 /// schemes. If these arguments are zero, simple binning schemes are
209 /// constructed which correspond to the axes of the histogram
210 /// <b>hist_A</b>.
211 ///
212 /// The parameters
213 /// <b>regmode, densityMode, regularisationDistribution, regularisationAxisSteering</b>
214 /// together control how the initial matrix L of regularisation conditions
215 /// is constructed. as explained in RegularizeDistribution().
216 
218 (const TH2 *hist_A, EHistMap histmap,ERegMode regmode,EConstraint constraint,
219  EDensityMode densityMode,const TUnfoldBinning *outputBins,
220  const TUnfoldBinning *inputBins,const char *regularisationDistribution,
221  const char *regularisationAxisSteering) :
222  TUnfoldSys(hist_A,histmap,kRegModeNone,constraint)
223 {
224  fRegularisationConditions=0;
225  // set up binning schemes
226  fConstOutputBins = outputBins;
227  fOwnedOutputBins = 0;
228  TAxis const *genAxis,*detAxis;
229  if(histmap==kHistMapOutputHoriz) {
230  genAxis=hist_A->GetXaxis();
231  detAxis=hist_A->GetYaxis();
232  } else {
233  genAxis=hist_A->GetYaxis();
234  detAxis=hist_A->GetXaxis();
235  }
236  if(!fConstOutputBins) {
237  // underflow and overflow are included in the binning scheme
238  // They may be used on generator level
239  fOwnedOutputBins =
240  new TUnfoldBinning(*genAxis,1,1);
241  fConstOutputBins = fOwnedOutputBins;
242  }
243  // check whether binning scheme is valid
244  if(fConstOutputBins->GetParentNode()) {
245  Error("TUnfoldDensity",
246  "Invalid output binning scheme (node is not the root node)");
247  }
248  fConstInputBins = inputBins;
249  fOwnedInputBins = 0;
250  if(!fConstInputBins) {
251  // underflow and overflow are not included in the binning scheme
252  // They are still used to count events which have not been reconstructed
253  fOwnedInputBins =
254  new TUnfoldBinning(*detAxis,0,0);
255  fConstInputBins = fOwnedInputBins;
256  }
257  if(fConstInputBins->GetParentNode()) {
258  Error("TUnfoldDensity",
259  "Invalid input binning scheme (node is not the root node)");
260  }
261  // check whether binning scheme matches with the histogram
262  // in terms of total number of bins
263  Int_t nOut=genAxis->GetNbins();
264  Int_t nOutMappedT=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins(kTRUE));
265  Int_t nOutMappedF=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins
266  (fOwnedOutputBins));
267  if((nOutMappedT!= nOut)&&(nOutMappedF!=nOut)) {
268  Error("TUnfoldDensity",
269  "Output binning incompatible number of bins: axis %d binning scheme %d (%d)",
270  nOut,nOutMappedT,nOutMappedF);
271  }
272  // check whether binning scheme matches with the histogram
273  Int_t nInput=detAxis->GetNbins();
274  Int_t nInputMappedT=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins(kTRUE));
275  Int_t nInputMappedF=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins
276  (fOwnedInputBins));
277  if((nInputMappedT!= nInput)&&(nInputMappedF!= nInput)) {
278  Error("TUnfoldDensity",
279  "Input binning incompatible number of bins:axis %d binning scheme %d (%d) ",
280  nInput,nInputMappedT,nInputMappedF);
281  }
282 
283  // report detailed list of excluded bins
284  for (Int_t ix = 0; ix <= nOut+1; ix++) {
285  if(fHistToX[ix]<0) {
286  Info("TUnfold","*NOT* unfolding bin %s",(char const *)GetOutputBinName(ix));
287  }
288  }
289 
290  // set up the regularisation here
291  if(regmode !=kRegModeNone) {
292  RegularizeDistribution
293  (regmode,densityMode,regularisationDistribution,
294  regularisationAxisSteering);
295  }
296 }
297 
298 ////////////////////////////////////////////////////////////////////////////////
299 /// Get bin name of an output bin.
300 ///
301 /// \param[in] iBinX bin number
302 ///
303 /// Return value: name of the bin. The name is constructed from the
304 /// entries in the binning scheme and includes information about the
305 /// bin borders etc.
306 
308  if(!fConstOutputBins) return TUnfold::GetOutputBinName(iBinX);
309  else return fConstOutputBins->GetBinName(iBinX);
310 }
311 
312 ////////////////////////////////////////////////////////////////////////////////
313 /// Density correction factor for a given bin.
314 ///
315 /// \param[in] densityMode type of factor to calculate
316 /// \param[in] iBin bin number
317 ///
318 /// return a multiplicative factor, for scaling the regularisation
319 /// conditions from this bin.
320 ///
321 /// For densityMode=kDensityModeNone the factor is set to unity.
322 /// For densityMode=kDensityModeBinWidth
323 /// the factor is set to 1/binArea
324 /// where the binArea is the product of the bin widths in all
325 /// dimensions. If the width of a bin is zero or can not be
326 /// determined, the factor is set to zero.
327 /// For densityMode=kDensityModeUser the factor is determined from the
328 /// method TUnfoldBinning::GetBinFactor().
329 /// For densityMode=kDensityModeBinWidthAndUser, the results of
330 /// kDensityModeBinWidth and kDensityModeUser are multiplied.
331 
333 (EDensityMode densityMode,Int_t iBin) const
334 {
335  Double_t factor=1.0;
336  if((densityMode == kDensityModeBinWidth)||
337  (densityMode == kDensityModeBinWidthAndUser)) {
338  Double_t binSize=fConstOutputBins->GetBinSize(iBin);
339  if(binSize>0.0) factor /= binSize;
340  else factor=0.0;
341  }
342  if((densityMode == kDensityModeUser)||
343  (densityMode == kDensityModeBinWidthAndUser)) {
344  factor *= fConstOutputBins->GetBinFactor(iBin);
345  }
346  return factor;
347 }
348 
349 ////////////////////////////////////////////////////////////////////////////////
350 /// Set up regularisation conditions.
351 ///
352 /// \param[in] regmode basic regularisation mode (see class TUnfold)
353 /// \param[in] densityMode how to apply bin-wise factors
354 /// \param[in] distribution name of the TUnfoldBinning node for which
355 /// the regularisation conditions shall be set (zero matches all nodes)
356 /// \param[in] axisSteering regularisation fine-tuning
357 ///
358 /// <b>axisSteering</b> is a string with several tokens, separated by
359 /// a semicolon: `"axisName[options];axisName[options];..."`.
360 ///
361 /// - <b>axisName</b>:
362 /// the name of an axis. The special name * matches all.
363 /// So the argument <b>distribution</b> selects one (or all)
364 /// distributions. Within the selected distribution(s), steering options may be
365 /// specified for each axis (or for all axes) to define the
366 /// regularisation conditions.
367 /// - <b>options</b>
368 /// one or several character as follows:
369 /// - u : exclude underflow bin from derivatives along this axis
370 /// - o : exclude overflow bin from derivatives along this axis
371 /// - U : exclude underflow bin
372 /// - O : exclude overflow bin
373 /// - b : use bin width for derivative calculation
374 /// - B : same as 'b', in addition normalize to average bin width
375 /// - N : completely exclude derivatives along this axis
376 /// - p : axis is periodic (e.g. azimuthal angle), so
377 /// include derivatives built from combinations involving bins at
378 /// both ends of the axis "wrap around"
379 ///
380 /// example: <b>axisSteering</b>=`"*[UOB]"` uses bin widths to calculate
381 /// derivatives but underflow/overflow bins are not regularized
382 
384 (ERegMode regmode,EDensityMode densityMode,const char *distribution,
385  const char *axisSteering)
386 {
387 
388  RegularizeDistributionRecursive(GetOutputBinning(),regmode,densityMode,
389  distribution,axisSteering);
390 }
391 
392 ////////////////////////////////////////////////////////////////////////////////
393 /// Recursively add regularisation conditions for this node and its children.
394 ///
395 /// \param[in] binning current node
396 /// \param[in] regmode regularisation mode
397 /// \param[in] densityMode type of regularisation scaling
398 /// \param[in] distribution target distribution(s) name
399 /// \param[in] axisSteering steering within the target distribution(s)
400 
402 (const TUnfoldBinning *binning,ERegMode regmode,
403  EDensityMode densityMode,const char *distribution,const char *axisSteering) {
404  if((!distribution)|| !TString(distribution).CompareTo(binning->GetName())) {
405  RegularizeOneDistribution(binning,regmode,densityMode,axisSteering);
406  }
407  for(const TUnfoldBinning *child=binning->GetChildNode();child;
408  child=child->GetNextNode()) {
409  RegularizeDistributionRecursive(child,regmode,densityMode,distribution,
410  axisSteering);
411  }
412 }
413 
414 ////////////////////////////////////////////////////////////////////////////////
415 /// Regularize the distribution of the given node.
416 ///
417 /// \param[in] binning current node
418 /// \param[in] regmode regularisation mode
419 /// \param[in] densityMode type of regularisation scaling
420 /// \param[in] axisSteering detailed steering for the axes of the distribution
421 
423 (const TUnfoldBinning *binning,ERegMode regmode,
424  EDensityMode densityMode,const char *axisSteering)
425 {
426 #ifdef DEBUG
427  cout<<"TUnfoldDensity::RegularizeOneDistribution node="
428  <<binning->GetName()<<" "<<regmode<<" "<<densityMode
429  <<" "<<(axisSteering ? axisSteering : "")<<"\n";
430 #endif
431  if(!fRegularisationConditions)
432  fRegularisationConditions=new TUnfoldBinning("regularisation");
433 
434  TUnfoldBinning *thisRegularisationBinning=
435  fRegularisationConditions->AddBinning(binning->GetName());
436 
437  // decode steering
438  Int_t isOptionGiven[8];
439  binning->DecodeAxisSteering(axisSteering,"uUoObBpN",isOptionGiven);
440  // U implies u
441  isOptionGiven[0] |= isOptionGiven[1];
442  // O implies o
443  isOptionGiven[2] |= isOptionGiven[3];
444  // B implies b
445  isOptionGiven[4] |= isOptionGiven[5];
446  // option N is removed if any other option is on
447  for(Int_t i=0;i<7;i++) {
448  isOptionGiven[7] &= ~isOptionGiven[i];
449  }
450  // option "c" does not work with options UuOo
451  if(isOptionGiven[6] & (isOptionGiven[0]|isOptionGiven[2]) ) {
452  Error("RegularizeOneDistribution",
453  "axis steering %s is not valid",axisSteering);
454  }
455 #ifdef DEBUG
456  cout<<" "<<isOptionGiven[0]
457  <<" "<<isOptionGiven[1]
458  <<" "<<isOptionGiven[2]
459  <<" "<<isOptionGiven[3]
460  <<" "<<isOptionGiven[4]
461  <<" "<<isOptionGiven[5]
462  <<" "<<isOptionGiven[6]
463  <<" "<<isOptionGiven[7]
464  <<"\n";
465 #endif
466  Info("RegularizeOneDistribution","regularizing %s regMode=%d"
467  " densityMode=%d axisSteering=%s",
468  binning->GetName(),(Int_t) regmode,(Int_t)densityMode,
469  axisSteering ? axisSteering : "");
470  Int_t startBin=binning->GetStartBin();
471  Int_t endBin=startBin+ binning->GetDistributionNumberOfBins();
472  std::vector<Double_t> factor(endBin-startBin);
473  Int_t nbin=0;
474  for(Int_t bin=startBin;bin<endBin;bin++) {
475  factor[bin-startBin]=GetDensityFactor(densityMode,bin);
476  if(factor[bin-startBin] !=0.0) nbin++;
477  }
478 #ifdef DEBUG
479  cout<<"initial number of bins "<<nbin<<"\n";
480 #endif
481  Int_t dimension=binning->GetDistributionDimension();
482 
483  // decide whether to skip underflow/overflow bins
484  nbin=0;
485  for(Int_t bin=startBin;bin<endBin;bin++) {
486  Int_t uStatus,oStatus;
487  binning->GetBinUnderflowOverflowStatus(bin,&uStatus,&oStatus);
488  if(uStatus & isOptionGiven[1]) factor[bin-startBin]=0.;
489  if(oStatus & isOptionGiven[3]) factor[bin-startBin]=0.;
490  if(factor[bin-startBin] !=0.0) nbin++;
491  }
492 #ifdef DEBUG
493  cout<<"after underflow/overflow bin removal "<<nbin<<"\n";
494 #endif
495  if(regmode==kRegModeSize) {
496  Int_t nRegBins=0;
497  // regularize all bins of the distribution, possibly excluding
498  // underflow/overflow bins
499  for(Int_t bin=startBin;bin<endBin;bin++) {
500  if(factor[bin-startBin]==0.0) continue;
501  if(AddRegularisationCondition(bin,factor[bin-startBin])) {
502  nRegBins++;
503  }
504  }
505  if(nRegBins) {
506  thisRegularisationBinning->AddBinning("size",nRegBins);
507  }
508  } else if((regmode==kRegModeDerivative)||(regmode==kRegModeCurvature)) {
509  for(Int_t direction=0;direction<dimension;direction++) {
510  // for each direction
511  Int_t nRegBins=0;
512  Int_t directionMask=(1<<direction);
513  if(isOptionGiven[7] & directionMask) {
514 #ifdef DEBUG
515  cout<<"skip direction "<<direction<<"\n";
516 #endif
517  continue;
518  }
519  Double_t binDistanceNormalisation=
520  (isOptionGiven[5] & directionMask) ?
522  (direction,isOptionGiven[0] & directionMask,
523  isOptionGiven[2] & directionMask) : 1.0;
524  for(Int_t bin=startBin;bin<endBin;bin++) {
525  // check whether bin is excluded
526  if(factor[bin-startBin]==0.0) continue;
527  // for each bin, find the neighbour bins
528  Int_t iPrev,iNext;
529  Double_t distPrev,distNext;
530  Int_t error=binning->GetBinNeighbours
531  (bin,direction,&iPrev,&distPrev,&iNext,&distNext,
532  isOptionGiven[6] & directionMask);
533  if(error) {
534  Error("RegularizeOneDistribution",
535  "invalid option %s (isPeriodic) for axis %s"
536  " (has underflow or overflow)",axisSteering,
537  binning->GetDistributionAxisLabel(direction).Data());
538  }
539  if((regmode==kRegModeDerivative)&&(iNext>=0)) {
540  Double_t f0 = -factor[bin-startBin];
541  Double_t f1 = factor[iNext-startBin];
542  if(isOptionGiven[4] & directionMask) {
543  if(distNext>0.0) {
544  f0 *= binDistanceNormalisation/distNext;
545  f1 *= binDistanceNormalisation/distNext;
546  } else {
547  f0=0.;
548  f1=0.;
549  }
550  }
551  if((f0==0.0)||(f1==0.0)) continue;
552  if(AddRegularisationCondition(bin,f0,iNext,f1)) {
553  nRegBins++;
554 #ifdef DEBUG
555  std::cout<<"Added Reg: bin "<<bin<<" "<<f0
556  <<" next: "<<iNext<<" "<<f1<<"\n";
557 #endif
558  }
559  } else if((regmode==kRegModeCurvature)&&(iPrev>=0)&&(iNext>=0)) {
560  Double_t f0 = factor[iPrev-startBin];
561  Double_t f1 = -factor[bin-startBin];
562  Double_t f2 = factor[iNext-startBin];
563  if(isOptionGiven[4] & directionMask) {
564  if((distPrev<0.)&&(distNext>0.)) {
565  distPrev= -distPrev;
566  Double_t f=TMath::Power(binDistanceNormalisation,2.)/
567  (distPrev+distNext);
568  f0 *= f/distPrev;
569  f1 *= f*(1./distPrev+1./distNext);
570  f2 *= f/distNext;
571  } else {
572  f0=0.;
573  f1=0.;
574  f2=0.;
575  }
576  }
577  if((f0==0.0)||(f1==0.0)||(f2==0.0)) continue;
578  if(AddRegularisationCondition(iPrev,f0,bin,f1,iNext,f2)) {
579  nRegBins++;
580 #ifdef DEBUG
581  std::cout<<"Added Reg: prev "<<iPrev<<" "<<f0
582  <<" bin: "<<bin<<" "<<f1
583  <<" next: "<<iNext<<" "<<f2<<"\n";
584 #endif
585  }
586  }
587  }
588  if(nRegBins) {
589  TString name;
590  if(regmode==kRegModeDerivative) {
591  name="derivative_";
592  } else if(regmode==kRegModeCurvature) {
593  name="curvature_";
594  }
595  name += binning->GetDistributionAxisLabel(direction);
596  thisRegularisationBinning->AddBinning(name,nRegBins);
597  }
598  }
599  }
600 #ifdef DEBUG
601  //fLsquared->Print();
602 #endif
603 }
604 
605 ///////////////////////////////////////////////////////////////////////
606 /// retrieve unfolding result as a new histogram
607 ///
608 /// \param[in] histogramName name of the histogram
609 /// \param[in] histogramTitle (default=0) title of the histogram
610 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
611 /// \param[in] axisSteering (default=0) detailed steering within the extracted
612 /// distribution
613 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
614 /// proper binning and axis labels
615 ///
616 /// return value: pointer to a new histogram. If
617 /// <b>useAxisBinning</b> is set and if the selected distribution fits
618 /// into a root histogram (1,2,3 dimensions) then return a histogram
619 /// with the proper binning on each axis. Otherwise, return a 1D
620 /// histogram with equidistant binning. If the histogram title is
621 /// zero, a title is assigned automatically.
622 ///
623 /// The <b>axisSteering</b> is defines as follows: "axis[mode];axis[mode];..."
624 /// where:
625 ///
626 /// - axis = name of an axis or *
627 /// - mode = any combination of the letters CUO0123456789
628 ///
629 /// - C collapse axis into one bin (add up results). If
630 /// any of the numbers 0-9 are given in addition, only these bins are added up.
631 /// Here bins are counted from zero, whereas in root, bins are counted
632 /// from 1. Obviously, this only works for up to 10 bins.
633 /// - U discard underflow bin
634 /// - O discard overflow bin
635 ///
636 /// examples: imagine the binning has two axis, named x and y.
637 ///
638 /// - "*[UO]" exclude underflow and overflow bins for all axis.
639 /// So here a TH2 is returned but all undeflow and overflow bins are empty
640 /// - "x[UOC123]" integrate over the variable x but only using the
641 /// bins 1,2,3 and not the underflow and overflow in x.
642 /// Here a TH1 is returned, the axis is labelled "y" and
643 /// the underflow and overflow (in y) are filled. However only the x-bins
644 /// 1,2,3 are used to determine the content.
645 /// - "x[C];y[UO]" integrate over the variable x, including
646 /// underflow and overflow but exclude underflow and overflow in y.
647 /// Here a TH1 is returned, the axis is labelled "y". The underflow
648 /// and overflow in y are empty.
649 
651 (const char *histogramName,const char *histogramTitle,
652  const char *distributionName,const char *axisSteering,
653  Bool_t useAxisBinning) const
654 {
655  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
656  Int_t *binMap=0;
657  TH1 *r=binning->CreateHistogram
658  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
659  if(r) {
660  TUnfoldSys::GetOutput(r,binMap);
661  }
662  if(binMap) {
663  delete [] binMap;
664  }
665  return r;
666 }
667 
668 ////////////////////////////////////////////////////////////////////////////////
669 /// Retrieve bias vector as a new histogram.
670 ///
671 /// \param[in] histogramName name of the histogram
672 /// \param[in] histogramTitle (default=0) title of the histogram
673 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
674 /// \param[in] axisSteering (default=0) detailed steering within the extracted
675 /// distribution
676 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
677 /// proper binning and axis labels
678 ///
679 /// returns a new histogram. See method GetOutput() for a detailed
680 /// description of the arguments
681 
683 (const char *histogramName,const char *histogramTitle,
684  const char *distributionName,const char *axisSteering,
685  Bool_t useAxisBinning) const
686 {
687  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
688  Int_t *binMap=0;
689  TH1 *r=binning->CreateHistogram
690  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
691  if(r) {
692  TUnfoldSys::GetBias(r,binMap);
693  }
694  if(binMap) delete [] binMap;
695  return r;
696 }
697 
698 ////////////////////////////////////////////////////////////////////////////////
699 /// Retrieve unfolding result folded back as a new histogram.
700 ///
701 /// \param[in] histogramName name of the histogram
702 /// \param[in] histogramTitle (default=0) title of the histogram
703 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
704 /// \param[in] axisSteering (default=0) detailed steering within the extracted
705 /// distribution
706 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
707 /// proper binning and axis labels
708 /// \param[in] addBgr (default=false) if true, include the background
709 /// contribution (useful for direct comparison to data)
710 ///
711 /// returns a new histogram. See method GetOutput() for a detailed
712 /// description of the arguments
713 
715 (const char *histogramName,const char *histogramTitle,
716  const char *distributionName,const char *axisSteering,
717  Bool_t useAxisBinning,Bool_t addBgr) const
718 {
719  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
720  Int_t *binMap=0;
721  TH1 *r=binning->CreateHistogram
722  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
723  if(r) {
724  TUnfoldSys::GetFoldedOutput(r,binMap);
725  if(addBgr) {
726  TUnfoldSys::GetBackground(r,0,binMap,0,kFALSE);
727  }
728  }
729  if(binMap) delete [] binMap;
730  return r;
731 }
732 
733 ////////////////////////////////////////////////////////////////////////////////
734 /// Retrieve a background source in a new histogram.
735 ///
736 /// \param[in] histogramName name of the histogram
737 /// \param[in] bgrSource the background source to retrieve
738 /// \param[in] histogramTitle (default=0) title of the histogram
739 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
740 /// \param[in] axisSteering (default=0) detailed steering within the extracted
741 /// distribution
742 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
743 /// proper binning and axis labels
744 /// \param[in] includeError (default=3) type of background errors to
745 /// be included (+1 uncorrelated bgr errors, +2 correlated bgr errors)
746 ///
747 /// returns a new histogram. See method GetOutput() for a detailed
748 /// description of the arguments
749 
751 (const char *histogramName,const char *bgrSource,const char *histogramTitle,
752  const char *distributionName,const char *axisSteering,Bool_t useAxisBinning,
753  Int_t includeError) const
754 {
755  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
756  Int_t *binMap=0;
757  TH1 *r=binning->CreateHistogram
758  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
759  if(r) {
760  TUnfoldSys::GetBackground(r,bgrSource,binMap,includeError,kTRUE);
761  }
762  if(binMap) delete [] binMap;
763  return r;
764 }
765 
766 ////////////////////////////////////////////////////////////////////////////////
767 /// Retrieve input distribution in a new histogram.
768 ///
769 /// \param[in] histogramName name of the histogram
770 /// \param[in] histogramTitle (default=0) title of the histogram
771 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
772 /// \param[in] axisSteering (default=0) detailed steering within the extracted
773 /// distribution
774 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
775 /// proper binning and axis labels
776 ///
777 /// returns a new histogram. See method GetOutput() for a detailed
778 /// description of the arguments
779 
781 (const char *histogramName,const char *histogramTitle,
782  const char *distributionName,const char *axisSteering,
783  Bool_t useAxisBinning) const
784 {
785  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
786  Int_t *binMap=0;
787  TH1 *r=binning->CreateHistogram
788  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
789  if(r) {
790  TUnfoldSys::GetInput(r,binMap);
791  }
792  if(binMap) delete [] binMap;
793  return r;
794 }
795 
796 ////////////////////////////////////////////////////////////////////////////////
797 /// Retrieve global correlation coefficients including all uncertainty sources.
798 ///
799 /// \param[in] histogramName name of the histogram
800 /// \param[in] histogramTitle (default=0) title of the histogram
801 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
802 /// \param[in] axisSteering (default=0) detailed steering within the extracted
803 /// distribution
804 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
805 /// proper binning and axis labels
806 /// \param[out] ematInv (default=0) to return the inverse covariance matrix
807 ///
808 /// returns a new histogram. See method GetOutput() for a detailed
809 /// description of the arguments. The inverse of the covariance matrix
810 /// is stored in a new histogram returned by <b>ematInv</b> if that
811 /// pointer is non-zero.
812 
814 (const char *histogramName,const char *histogramTitle,
815  const char *distributionName,const char *axisSteering,
816  Bool_t useAxisBinning,TH2 **ematInv) {
817  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
818  Int_t *binMap=0;
819  TH1 *r=binning->CreateHistogram
820  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
821  if(r) {
822  TH2 *invEmat=0;
823  if(ematInv) {
824  if(r->GetDimension()==1) {
825  TString ematName(histogramName);
826  ematName += "_inverseEMAT";
827  Int_t *binMap2D=0;
828  invEmat=binning->CreateErrorMatrixHistogram
829  (ematName,useAxisBinning,&binMap2D,histogramTitle,
830  axisSteering);
831  if(binMap2D) delete [] binMap2D;
832  } else {
833  Error("GetRhoItotal",
834  "can not return inverse of error matrix for this binning");
835  }
836  }
837  TUnfoldSys::GetRhoItotal(r,binMap,invEmat);
838  if(invEmat) {
839  *ematInv=invEmat;
840  }
841  }
842  if(binMap) delete [] binMap;
843  return r;
844 }
845 
846 ////////////////////////////////////////////////////////////////////////////////
847 /// Retrieve global correlation coefficients including input
848 /// (statistical) and background uncertainties.
849 ///
850 /// \param[in] histogramName name of the histogram
851 /// \param[in] histogramTitle (default=0) title of the histogram
852 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
853 /// \param[in] axisSteering (default=0) detailed steering within the extracted
854 /// distribution
855 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
856 /// proper binning and axis labels
857 /// \param[out] ematInv (default=0) to return the inverse covariance matrix
858 ///
859 /// returns a new histogram. See method GetOutput() for a detailed
860 /// description of the arguments. The inverse of the covariance matrix
861 /// is stored in a new histogram returned by <b>ematInv</b> if that
862 /// pointer is non-zero.
863 
865 (const char *histogramName,const char *histogramTitle,
866  const char *distributionName,const char *axisSteering,
867  Bool_t useAxisBinning,TH2 **ematInv) {
868  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
869  Int_t *binMap=0;
870  TH1 *r=binning->CreateHistogram
871  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
872  if(r) {
873  TH2 *invEmat=0;
874  if(ematInv) {
875  if(r->GetDimension()==1) {
876  TString ematName(histogramName);
877  ematName += "_inverseEMAT";
878  Int_t *binMap2D=0;
879  invEmat=binning->CreateErrorMatrixHistogram
880  (ematName,useAxisBinning,&binMap2D,histogramTitle,
881  axisSteering);
882  if(binMap2D) delete [] binMap2D;
883  } else {
884  Error("GetRhoItotal",
885  "can not return inverse of error matrix for this binning");
886  }
887  }
888  TUnfoldSys::GetRhoI(r,binMap,invEmat);
889  if(invEmat) {
890  *ematInv=invEmat;
891  }
892  }
893  if(binMap) delete [] binMap;
894  return r;
895 }
896 
897 ////////////////////////////////////////////////////////////////////////////////
898 /// Retrieve a correlated systematic 1-sigma shift.
899 ///
900 /// \param[in] source identifier of the systematic uncertainty source
901 /// \param[in] histogramName name of the histogram
902 /// \param[in] histogramTitle (default=0) title of the histogram
903 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
904 /// \param[in] axisSteering (default=0) detailed steering within the extracted
905 /// distribution
906 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
907 /// proper binning and axis labels
908 ///
909 /// returns a new histogram. See method GetOutput() for a detailed
910 /// description of the arguments
911 
913 (const char *source,const char *histogramName,
914  const char *histogramTitle,const char *distributionName,
915  const char *axisSteering,Bool_t useAxisBinning) {
916  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
917  Int_t *binMap=0;
918  TH1 *r=binning->CreateHistogram
919  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
920  if(r) {
921  if(!TUnfoldSys::GetDeltaSysSource(r,source,binMap)) {
922  delete r;
923  r=0;
924  }
925  }
926  if(binMap) delete [] binMap;
927  return r;
928 }
929 
930 ////////////////////////////////////////////////////////////////////////////////
931 /// Retrieve systematic 1-sigma shift corresponding to a background
932 /// scale uncertainty.
933 ///
934 /// \param[in] bgrSource identifier of the background
935 /// \param[in] histogramName name of the histogram
936 /// \param[in] histogramTitle (default=0) title of the histogram
937 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
938 /// \param[in] axisSteering (default=0) detailed steering within the extracted
939 /// distribution
940 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
941 /// proper binning and axis labels
942 ///
943 /// returns a new histogram. See method GetOutput() for a detailed
944 /// description of the arguments
945 
947 (const char *bgrSource,const char *histogramName,
948  const char *histogramTitle,const char *distributionName,
949  const char *axisSteering,Bool_t useAxisBinning) {
950  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
951  Int_t *binMap=0;
952  TH1 *r=binning->CreateHistogram
953  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
954  if(r) {
955  if(!TUnfoldSys::GetDeltaSysBackgroundScale(r,bgrSource,binMap)) {
956  delete r;
957  r=0;
958  }
959  }
960  if(binMap) delete [] binMap;
961  return r;
962 }
963 
964 ////////////////////////////////////////////////////////////////////////////////
965 /// Retrieve 1-sigma shift corresponding to the previously specified uncertainty
966 /// on tau.
967 ///
968 /// \param[in] histogramName name of the histogram
969 /// \param[in] histogramTitle (default=0) title of the histogram
970 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
971 /// \param[in] axisSteering (default=0) detailed steering within the extracted
972 /// distribution
973 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
974 /// proper binning and axis labels
975 ///
976 /// returns a new histogram. See method GetOutput() for a detailed
977 /// description of the arguments
978 
980 (const char *histogramName,const char *histogramTitle,
981  const char *distributionName,const char *axisSteering,Bool_t useAxisBinning)
982 {
983  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
984  Int_t *binMap=0;
985  TH1 *r=binning->CreateHistogram
986  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
987  if(r) {
988  if(!TUnfoldSys::GetDeltaSysTau(r,binMap)) {
989  delete r;
990  r=0;
991  }
992  }
993  if(binMap) delete [] binMap;
994  return r;
995 }
996 
997 ////////////////////////////////////////////////////////////////////////////////
998 /// Retrieve correlation coefficients, including all uncertainties.
999 ///
1000 /// \param[in] histogramName name of the histogram
1001 /// \param[in] histogramTitle (default=0) title of the histogram
1002 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1003 /// \param[in] axisSteering (default=0) detailed steering within the extracted
1004 /// distribution
1005 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1006 /// proper binning and axis labels
1007 ///
1008 /// returns a new histogram. See method GetOutput() for a detailed
1009 /// description of the arguments
1010 
1012 (const char *histogramName,const char *histogramTitle,
1013  const char *distributionName,const char *axisSteering,
1014  Bool_t useAxisBinning)
1015 {
1016  TH2 *r=GetEmatrixTotal
1017  (histogramName,histogramTitle,distributionName,
1018  axisSteering,useAxisBinning);
1019  if(r) {
1020  for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
1021  Double_t e_i=r->GetBinContent(i,i);
1022  if(e_i>0.0) e_i=TMath::Sqrt(e_i);
1023  else e_i=0.0;
1024  for(Int_t j=0;j<=r->GetNbinsY()+1;j++) {
1025  if(i==j) continue;
1026  Double_t e_j=r->GetBinContent(j,j);
1027  if(e_j>0.0) e_j=TMath::Sqrt(e_j);
1028  else e_j=0.0;
1029  Double_t e_ij=r->GetBinContent(i,j);
1030  if((e_i>0.0)&&(e_j>0.0)) {
1031  r->SetBinContent(i,j,e_ij/e_i/e_j);
1032  } else {
1033  r->SetBinContent(i,j,0.0);
1034  }
1035  }
1036  }
1037  for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
1038  if(r->GetBinContent(i,i)>0.0) {
1039  r->SetBinContent(i,i,1.0);
1040  } else {
1041  r->SetBinContent(i,i,0.0);
1042  }
1043  }
1044  }
1045  return r;
1046 }
1047 
1048 ////////////////////////////////////////////////////////////////////////////////
1049 /// Retrieve covariance contribution from uncorrelated (statistical)
1050 /// uncertainties of the response matrix.
1051 ///
1052 /// \param[in] histogramName name of the histogram
1053 /// \param[in] histogramTitle (default=0) title of the histogram
1054 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1055 /// \param[in] axisSteering (default=0) detailed steering within the extracted
1056 /// distribution
1057 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1058 /// proper binning and axis labels
1059 ///
1060 /// returns a new histogram. See method GetOutput() for a detailed
1061 /// description of the arguments
1062 
1064 (const char *histogramName,const char *histogramTitle,
1065  const char *distributionName,const char *axisSteering,
1066  Bool_t useAxisBinning)
1067 {
1068  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1069  Int_t *binMap=0;
1070  TH2 *r=binning->CreateErrorMatrixHistogram
1071  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1072  if(r) {
1074  }
1075  if(binMap) delete [] binMap;
1076  return r;
1077 }
1078 
1079 ////////////////////////////////////////////////////////////////////////////////
1080 /// Retrieve covariance contribution from uncorrelated background uncertainties.
1081 ///
1082 /// \param[in] bgrSource identifier of the background
1083 /// \param[in] histogramName name of the histogram
1084 /// \param[in] histogramTitle (default=0) title of the histogram
1085 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1086 /// \param[in] axisSteering (default=0) detailed steering within the extracted
1087 /// distribution
1088 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1089 /// proper binning and axis labels
1090 ///
1091 /// returns a new histogram. See method GetOutput() for a detailed
1092 /// description of the arguments
1093 
1095 (const char *bgrSource,const char *histogramName,
1096  const char *histogramTitle,const char *distributionName,
1097  const char *axisSteering,Bool_t useAxisBinning)
1098 {
1099  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1100  Int_t *binMap=0;
1101  TH2 *r=binning->CreateErrorMatrixHistogram
1102  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1103  if(r) {
1105  }
1106  if(binMap) delete [] binMap;
1107  return r;
1108 }
1109 
1110 ////////////////////////////////////////////////////////////////////////////////
1111 /// Get covariance contribution from the input uncertainties (data
1112 /// statistical uncertainties).
1113 ///
1114 /// \param[in] histogramName name of the histogram
1115 /// \param[in] histogramTitle (default=0) title of the histogram
1116 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1117 /// \param[in] axisSteering (default=0) detailed steering within the extracted
1118 /// distribution
1119 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1120 /// proper binning and axis labels
1121 ///
1122 /// returns a new histogram. See method GetOutput() for a detailed
1123 /// description of the arguments
1124 
1126 (const char *histogramName,const char *histogramTitle,
1127  const char *distributionName,const char *axisSteering,
1128  Bool_t useAxisBinning)
1129 {
1130  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1131  Int_t *binMap=0;
1132  TH2 *r=binning->CreateErrorMatrixHistogram
1133  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1134  if(r) {
1135  TUnfoldSys::GetEmatrixInput(r,binMap);
1136  }
1137  if(binMap) delete [] binMap;
1138  return r;
1139 }
1140 
1141 ////////////////////////////////////////////////////////////////////////////////
1142 /// Get matrix of probabilities in a new histogram.
1143 ///
1144 /// \param[in] histogramName name of the histogram
1145 /// \param[in] histogramTitle (default=0) title of the histogram
1146 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1147 /// proper binning and axis labels
1148 ///
1149 /// returns a new histogram. if histogramTitle is null, choose a title
1150 /// automatically.
1151 
1153 (const char *histogramName,const char *histogramTitle,
1154  Bool_t useAxisBinning) const
1155 {
1157  (fConstOutputBins,fConstInputBins,histogramName,
1158  useAxisBinning,useAxisBinning,histogramTitle);
1159  TUnfold::GetProbabilityMatrix(r,kHistMapOutputHoriz);
1160  return r;
1161 }
1162 
1163 ////////////////////////////////////////////////////////////////////////////////
1164 /// Get covariance matrix including all contributions.
1165 ///
1166 /// \param[in] histogramName name of the histogram
1167 /// \param[in] histogramTitle (default=0) title of the histogram
1168 /// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1169 /// \param[in] axisSteering (default=0) detailed steering within the extracted
1170 /// distribution
1171 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1172 /// proper binning and axis labels
1173 ///
1174 /// returns a new histogram. See method GetOutput() for a detailed
1175 /// description of the arguments
1176 
1178 (const char *histogramName,const char *histogramTitle,
1179  const char *distributionName,const char *axisSteering,
1180  Bool_t useAxisBinning)
1181 {
1182  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1183  Int_t *binMap=0;
1184  TH2 *r=binning->CreateErrorMatrixHistogram
1185  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1186  if(r) {
1187  TUnfoldSys::GetEmatrixTotal(r,binMap);
1188  }
1189  if(binMap) delete [] binMap;
1190  return r;
1191 }
1192 
1193 ////////////////////////////////////////////////////////////////////////////////
1194 /// Access matrix of regularisation conditions in a new histogram.
1195 ///
1196 /// \param[in] histogramName name of the histogram
1197 /// \param[in] histogramTitle (default=0) title of the histogram
1198 /// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1199 /// proper binning and axis labels
1200 ///
1201 /// returns a new histogram. if histogramTitle is null, choose a title
1202 /// automatically.
1203 
1205 (const char *histogramName,const char *histogramTitle,Bool_t useAxisBinning)
1206 {
1207  if(fRegularisationConditions &&
1208  (fRegularisationConditions->GetEndBin()-
1209  fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
1210  Warning("GetL",
1211  "remove invalid scheme of regularisation conditions %d %d",
1212  fRegularisationConditions->GetEndBin(),fL->GetNrows());
1213  delete fRegularisationConditions;
1214  fRegularisationConditions=0;
1215  }
1216  if(!fRegularisationConditions) {
1217  fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1218  Warning("GetL","create flat regularisation conditions scheme");
1219  }
1221  (fConstOutputBins,fRegularisationConditions,histogramName,
1222  useAxisBinning,useAxisBinning,histogramTitle);
1223  TUnfold::GetL(r);
1224  return r;
1225 }
1226 
1227 ////////////////////////////////////////////////////////////////////////////////
1228 /// Get regularisation conditions multiplied by result vector minus bias
1229 /// L(x-biasScale*biasVector).
1230 ///
1231 /// \param[in] histogramName name of the histogram
1232 /// \param[in] histogramTitle (default=0) title of the histogram
1233 ///
1234 /// returns a new histogram.
1235 /// This is a measure of the level of regularization required per
1236 /// regularisation condition.
1237 /// If there are (negative or positive) spikes,
1238 /// these regularisation conditions dominate
1239 /// over the other regularisation conditions and may introduce
1240 /// the largest biases.
1241 
1243 (const char *histogramName,const char *histogramTitle)
1244 {
1245  TMatrixD dx(*GetX(), TMatrixD::kMinus, fBiasScale * (*fX0));
1246  TMatrixDSparse *Ldx=MultiplyMSparseM(fL,&dx);
1247  if(fRegularisationConditions &&
1248  (fRegularisationConditions->GetEndBin()-
1249  fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
1250  Warning("GetLxMinusBias",
1251  "remove invalid scheme of regularisation conditions %d %d",
1252  fRegularisationConditions->GetEndBin(),fL->GetNrows());
1253  delete fRegularisationConditions;
1254  fRegularisationConditions=0;
1255  }
1256  if(!fRegularisationConditions) {
1257  fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1258  Warning("GetLxMinusBias","create flat regularisation conditions scheme");
1259  }
1260  TH1 *r=fRegularisationConditions->CreateHistogram
1261  (histogramName,kFALSE,0,histogramTitle);
1262  const Int_t *Ldx_rows=Ldx->GetRowIndexArray();
1263  const Double_t *Ldx_data=Ldx->GetMatrixArray();
1264  for(Int_t row=0;row<Ldx->GetNrows();row++) {
1265  if(Ldx_rows[row]<Ldx_rows[row+1]) {
1266  r->SetBinContent(row+1,Ldx_data[Ldx_rows[row]]);
1267  }
1268  }
1269  delete Ldx;
1270  return r;
1271 }
1272 
1273 ////////////////////////////////////////////////////////////////////////////////
1274 /// Locate a binning node for the input (measured) quantities.
1275 ///
1276 /// \param[in] distributionName (default=0) distribution to look
1277 /// for. if zero, return the root node
1278 ///
1279 /// returns: pointer to a TUnfoldBinning object or zero if not found
1280 
1282 (const char *distributionName) const
1283 {
1284  // find binning scheme, input bins
1285  // distributionName : the distribution to locate
1286  return fConstInputBins->FindNode(distributionName);
1287 }
1288 
1289 ////////////////////////////////////////////////////////////////////////////////
1290 /// Locate a binning node for the unfolded (truth level) quantities.
1291 ///
1292 /// \param[in] distributionName (default=0) distribution to look
1293 /// for. if zero, return the root node
1294 ///
1295 /// returns: pointer to a TUnfoldBinning object or zero if not found
1296 
1298 (const char *distributionName) const
1299 {
1300  // find binning scheme, output bins
1301  // distributionName : the distribution to locate
1302  return fConstOutputBins->FindNode(distributionName);
1303 }
1304 
1305 ////////////////////////////////////////////////////////////////////////////////
1306 /// Scan a function wrt tau and determine the minimum.
1307 ///
1308 /// \param[in] nPoint number of points to be scanned
1309 /// \param[in] tauMin smallest tau value to study
1310 /// \param[in] tauMax largest tau value to study
1311 /// \param[out] scanResult the scanned function wrt log(tau)
1312 /// \param[in] mode 1st parameter for the scan function
1313 /// \param[in] distribution 2nd parameter for the scan function
1314 /// \param[in] projectionMode 3rd parameter for the scan function
1315 /// \param[out] lCurvePlot for monitoring, shows the L-curve
1316 /// \param[out] logTauXPlot for monitoring, L-curve(X) as a function of log(tau)
1317 /// \param[out] logTauYPlot for monitoring, L-curve(Y) as a function of log(tau)
1318 ///
1319 /// Return value: the coordinate number on the curve <b>scanResult</b>
1320 /// which corresponds to the minimum
1321 ///
1322 /// The function is scanned by repeating the following steps <b>nPoint</b>
1323 /// times
1324 ///
1325 /// 1. Choose a value of tau
1326 /// 2. Perform the unfolding for this choice of tau, DoUnfold(tau)
1327 /// 3. Determine the scan variable GetScanVariable()
1328 ///
1329 /// The method GetScanVariable() defines scans of correlation
1330 /// coefficients, where <b>mode</b> is chosen from the enum
1331 /// EScanTauMode. In addition one may set <b>distribution</b>
1332 /// and/or <b>projectionMode</b> to refine the calculation of
1333 /// correlations (e.g. restrict the calculation to the signal
1334 /// distribution and/or exclude underflow and overflow bins).
1335 /// See the documentation of GetScanVariable() for details.
1336 /// Alternative scan variables may be defined by overriding the
1337 /// GetScanVariable() method.
1338 ///
1339 /// Automatic choice of scan range: if (tauMin,tauMax) do not
1340 /// correspond to a valid tau range (e.g. tauMin=tauMax=0.0) then
1341 /// the tau range is determined automatically. Use with care!
1342 
1344 (Int_t nPoint,Double_t tauMin,Double_t tauMax,TSpline **scanResult,
1345  Int_t mode,const char *distribution,const char *axisSteering,
1346  TGraph **lCurvePlot,TSpline **logTauXPlot,TSpline **logTauYPlot)
1347 {
1348  typedef std::map<Double_t,Double_t> TauScan_t;
1349  typedef std::map<Double_t,std::pair<Double_t,Double_t> > LCurve_t;
1350  TauScan_t curve;
1351  LCurve_t lcurve;
1352 
1353  //==========================================================
1354  // algorithm:
1355  // (1) do the unfolding for nPoint-1 points
1356  // and store the results in the map
1357  // curve
1358  // (1a) store minimum and maximum tau to curve
1359  // (1b) insert additional points, until nPoint-1 values
1360  // have been calculated
1361  //
1362  // (2) determine the best choice of tau
1363  // do the unfolding for this point
1364  // and store the result in
1365  // curve
1366  // (3) return the result in scanResult
1367 
1368  //==========================================================
1369  // (1) do the unfolding for nPoint-1 points
1370  // and store the results in
1371  // curve
1372  // (1a) store minimum and maximum tau to curve
1373 
1374  if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
1375  // here no range is given, has to be determined automatically
1376  // the maximum tau is determined from the chi**2 values
1377  // observed from unfolding without regularization
1378 
1379  // first unfolding, without regularisation
1380  DoUnfold(0.0);
1381 
1382  // if the number of degrees of freedom is too small, create an error
1383  if(GetNdf()<=0) {
1384  Error("ScanTau","too few input bins, NDF<=0 %d",GetNdf());
1385  }
1386  Double_t X0=GetLcurveX();
1387  Double_t Y0=GetLcurveY();
1388  Double_t y0=GetScanVariable(mode,distribution,axisSteering);
1389  Info("ScanTau","logtau=-Infinity y=%lf X=%lf Y=%lf",y0,X0,Y0);
1390  {
1391  // unfolding guess maximum tau and store it
1392  Double_t logTau=
1393  0.5*(TMath::Log10(GetChi2A()+3.*TMath::Sqrt(GetNdf()+1.0))
1394  -GetLcurveY());
1395  DoUnfold(TMath::Power(10.,logTau));
1396  if((!TMath::Finite(GetLcurveX())) ||(!TMath::Finite(GetLcurveY()))) {
1397  Fatal("ScanTau","problem (missing regularisation?) X=%f Y=%f",
1398  GetLcurveX(),GetLcurveY());
1399  }
1400  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1401  curve[logTau]=y;
1402  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1403  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1404  GetLcurveX(),GetLcurveY());
1405  }
1406  // minimum tau is chosen such that it is less than
1407  // 1% different from the case of no regularisation
1408  // here, several points are inserted as needed
1409  while(((int)curve.size()<nPoint-1)&&
1410  ((TMath::Abs(GetLcurveX()-X0)>0.00432)||
1411  (TMath::Abs(GetLcurveY()-Y0)>0.00432))) {
1412  Double_t logTau=(*curve.begin()).first-0.5;
1413  DoUnfold(TMath::Power(10.,logTau));
1414  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1415  curve[logTau]=y;
1416  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1417  Info("ScanTay","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1418  GetLcurveX(),GetLcurveY());
1419  }
1420  } else {
1421  Double_t logTauMin=TMath::Log10(tauMin);
1422  Double_t logTauMax=TMath::Log10(tauMax);
1423  if(nPoint>1) {
1424  // insert maximum tau
1425  DoUnfold(TMath::Power(10.,logTauMax));
1426  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1427  curve[logTauMax]=y;
1428  lcurve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
1429  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMax,y,
1430  GetLcurveX(),GetLcurveY());
1431  }
1432  // insert minimum tau
1433  DoUnfold(TMath::Power(10.,logTauMin));
1434  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1435  curve[logTauMin]=y;
1436  lcurve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
1437  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMin,y,
1438  GetLcurveX(),GetLcurveY());
1439  }
1440 
1441  //==========================================================
1442  // (1b) insert additional points, until nPoint-1 values
1443  // have been calculated
1444  while((int)curve.size()<nPoint-1) {
1445  // insert additional points
1446  // points are inserted such that the largest interval in log(tau)
1447  // is divided into two smaller intervals
1448  // however, there is a penalty term for subdividing intervals
1449  // which are far away from the minimum
1450  TauScan_t::const_iterator i0,i1;
1451  i0=curve.begin();
1452  // locate minimum
1453  Double_t logTauYMin=(*i0).first;
1454  Double_t yMin=(*i0).second;
1455  for (; i0 != curve.end(); ++i0) {
1456  if((*i0).second<yMin) {
1457  yMin=(*i0).second;
1458  logTauYMin=(*i0).first;
1459  }
1460  }
1461  // insert additional points such that large log(tau) intervals
1462  // near the minimum rho are divided into two
1463  i0=curve.begin();
1464  i1=i0;
1465  Double_t distMax=0.0;
1466  Double_t logTau=0.0;
1467  for (++i1; i1 != curve.end(); ++i1) {
1468  Double_t dist;
1469  // check size of rho interval
1470  dist=TMath::Abs((*i0).first-(*i1).first)
1471  // penalty term if distance from rhoMax is large
1472  +0.25*TMath::Power(0.5*((*i0).first+(*i1).first)-logTauYMin,2.)/
1473  ((*curve.rbegin()).first-(*curve.begin()).first)/nPoint;
1474  if((dist<=0.0)||(dist>distMax)) {
1475  distMax=dist;
1476  logTau=0.5*((*i0).first+(*i1).first);
1477  }
1478  i0=i1;
1479  }
1480  DoUnfold(TMath::Power(10.,logTau));
1481  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1482  curve[logTau]=y;
1483  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1484  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1485  GetLcurveX(),GetLcurveY());
1486  }
1487 
1488  //==========================================================
1489  // (2) determine the best choice of tau
1490  // do the unfolding for this point
1491  // and store the result in
1492  // curve
1493 
1494  Double_t cTmin=0.0;
1495  {
1496  Double_t *cTi=new Double_t[curve.size()];
1497  Double_t *cCi=new Double_t[curve.size()];
1498  Int_t n=0;
1499  for (TauScan_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
1500  cTi[n]=(*i).first;
1501  cCi[n]=(*i).second;
1502  n++;
1503  }
1504  // create rho Spline
1505  TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n);
1506  // find the maximum of the curvature
1507  // if the parameter iskip is non-zero, then iskip points are
1508  // ignored when looking for the largest curvature
1509  // (there are problems with the curvature determined from the first
1510  // few points of splineX,splineY in the algorithm above)
1511  Int_t iskip=0;
1512  if(n>3) iskip=1;
1513  if(n>6) iskip=2;
1514  Double_t cCmin=cCi[iskip];
1515  cTmin=cTi[iskip];
1516  for(Int_t i=iskip;i<n-1-iskip;i++) {
1517  // find minimum on this spline section
1518  // check boundary conditions for x[i+1]
1519  Double_t xMin=cTi[i+1];
1520  Double_t yMin=cCi[i+1];
1521  if(cCi[i]<yMin) {
1522  yMin=cCi[i];
1523  xMin=cTi[i];
1524  }
1525  // find minimum for x[i]<x<x[i+1]
1526  // get spline coefficients and solve equation
1527  // derivative(x)==0
1528  Double_t x,y,b,c,d;
1529  splineC->GetCoeff(i,x,y,b,c,d);
1530  // coefficients of quadratic equation
1531  Double_t m_p_half=-c/(3.*d);
1532  Double_t q=b/(3.*d);
1533  Double_t discr=m_p_half*m_p_half-q;
1534  if(discr>=0.0) {
1535  // solution found
1536  discr=TMath::Sqrt(discr);
1537  Double_t xx;
1538  if(m_p_half>0.0) {
1539  xx = m_p_half + discr;
1540  } else {
1541  xx = m_p_half - discr;
1542  }
1543  Double_t dx=cTi[i+1]-x;
1544  // check first solution
1545  if((xx>0.0)&&(xx<dx)) {
1546  y=splineC->Eval(x+xx);
1547  if(y<yMin) {
1548  yMin=y;
1549  xMin=x+xx;
1550  }
1551  }
1552  // second solution
1553  if(xx !=0.0) {
1554  xx= q/xx;
1555  } else {
1556  xx=0.0;
1557  }
1558  // check second solution
1559  if((xx>0.0)&&(xx<dx)) {
1560  y=splineC->Eval(x+xx);
1561  if(y<yMin) {
1562  yMin=y;
1563  xMin=x+xx;
1564  }
1565  }
1566  }
1567  // check whether this local minimum is a global minimum
1568  if(yMin<cCmin) {
1569  cCmin=yMin;
1570  cTmin=xMin;
1571  }
1572  }
1573  delete splineC;
1574  delete[] cTi;
1575  delete[] cCi;
1576  }
1577  Double_t logTauFin=cTmin;
1578  DoUnfold(TMath::Power(10.,logTauFin));
1579  {
1580  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1581  curve[logTauFin]=y;
1582  lcurve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
1583  Info("ScanTau","Result logtau=%lf y=%lf X=%lf Y=%lf",logTauFin,y,
1584  GetLcurveX(),GetLcurveY());
1585  }
1586  //==========================================================
1587  // (3) return the result in
1588  // scanResult lCurve logTauX logTauY
1589 
1590  Int_t bestChoice=-1;
1591  if(curve.size()>0) {
1592  Double_t *y=new Double_t[curve.size()];
1593  Double_t *logT=new Double_t[curve.size()];
1594  int n=0;
1595  for (TauScan_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
1596  if(logTauFin==(*i).first) {
1597  bestChoice=n;
1598  }
1599  y[n]=(*i).second;
1600  logT[n]=(*i).first;
1601  n++;
1602  }
1603  if(scanResult) {
1604  TString name;
1605  name = TString::Format("scan(%d,",mode);
1606  if(distribution) name+= distribution;
1607  name += ",";
1608  if(axisSteering) name += axisSteering;
1609  name +=")";
1610  (*scanResult)=new TSpline3(name+"%log(tau)",logT,y,n);
1611  }
1612  delete[] y;
1613  delete[] logT;
1614  }
1615  if(lcurve.size()) {
1616  Double_t *logT=new Double_t[lcurve.size()];
1617  Double_t *x=new Double_t[lcurve.size()];
1618  Double_t *y=new Double_t[lcurve.size()];
1619  Int_t n=0;
1620  for (LCurve_t::const_iterator i = lcurve.begin(); i != lcurve.end(); ++i) {
1621  logT[n]=(*i).first;
1622  x[n]=(*i).second.first;
1623  y[n]=(*i).second.second;
1624  //cout<<logT[n]<<" "<< x[n]<<" "<<y[n]<<"\n";
1625  n++;
1626  }
1627  if(lCurvePlot) {
1628  *lCurvePlot=new TGraph(n,x,y);
1629  (*lCurvePlot)->SetTitle("L curve");
1630  }
1631  if(logTauXPlot)
1632  *logTauXPlot=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
1633  if(logTauYPlot)
1634  *logTauYPlot=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
1635  delete [] y;
1636  delete [] x;
1637  delete [] logT;
1638  }
1639  return bestChoice;
1640 }
1641 
1642 ////////////////////////////////////////////////////////////////////////////////
1643 /// Calculate the function for ScanTau().
1644 ///
1645 /// \param[in] mode the variable to be calculated
1646 /// \param[in] distribution distribution for which the variable
1647 /// is to be calculated
1648 /// \param[in] axisSteering detailed steering for selecting bins on
1649 /// the axes of the distribution (see method GetRhoItotal())
1650 ///
1651 /// return value: the scan result for the given choice of tau (for
1652 /// which the unfolding was performed prior to calling this method)
1653 ///
1654 /// In ScanTau() the unfolding is repeated for various choices of tau.
1655 /// For each tau, after unfolding, GetScanVariable() is called to
1656 /// determine the scan result for this choice of tau.
1657 ///
1658 /// the following modes are implemented
1659 ///
1660 /// - kEScanTauRhoAvg : average (stat+bgr) global correlation
1661 /// - kEScanTauRhoSquaredAvg : average (stat+bgr) global correlation squared
1662 /// - kEScanTauRhoMax : maximum (stat+bgr) global correlation
1663 /// - kEScanTauRhoAvgSys : average (stat+bgr+sys) global correlation
1664 /// - kEScanTauRhoAvgSquaredSys : average (stat+bgr+sys) global correlation squared
1665 /// - kEScanTauRhoMaxSys : maximum (stat+bgr+sys) global correlation
1666 
1668 (Int_t mode,const char *distribution,const char *axisSteering)
1669 {
1670  Double_t r=0.0;
1671  TString name("GetScanVariable(");
1672  name += TString::Format("%d,",mode);
1673  if(distribution) name += distribution;
1674  name += ",";
1675  if(axisSteering) name += axisSteering;
1676  name += ")";
1677  TH1 *rhoi=0;
1678  if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoMax)||
1679  (mode==kEScanTauRhoSquareAvg)) {
1680  rhoi=GetRhoIstatbgr(name,0,distribution,axisSteering,kFALSE);
1681  } else if((mode==kEScanTauRhoAvgSys)||(mode==kEScanTauRhoMaxSys)||
1682  (mode==kEScanTauRhoSquareAvgSys)) {
1683  rhoi=GetRhoItotal(name,0,distribution,axisSteering,kFALSE);
1684  }
1685  if(rhoi) {
1686  Double_t sum=0.0;
1687  Double_t sumSquare=0.0;
1688  Double_t rhoMax=0.0;
1689  Int_t n=0;
1690  for(Int_t i=0;i<=rhoi->GetNbinsX()+1;i++) {
1691  Double_t c=rhoi->GetBinContent(i);
1692  if(c>=0.) {
1693  if(c>rhoMax) rhoMax=c;
1694  sum += c;
1695  sumSquare += c*c;
1696  n ++;
1697  }
1698  }
1699  if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoAvgSys)) {
1700  r=sum/n;
1701  } else if((mode==kEScanTauRhoSquareAvg)||
1702  (mode==kEScanTauRhoSquareAvgSys)) {
1703  r=sum/n;
1704  } else {
1705  r=rhoMax;
1706  }
1707  // cout<<r<<" "<<GetRhoAvg()<<" "<<GetRhoMax()<<"\n";
1708  delete rhoi;
1709  } else {
1710  Fatal("GetScanVariable","mode %d not implemented",mode);
1711  }
1712  return r;
1713 }
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:47
virtual const Element * GetMatrixArray() const
TUnfoldBinning const * GetChildNode(void) const
first daughter node
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
static long int sum(long int i)
Definition: Factory.cxx:2258
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
Create a THxx histogram capable to hold the bins of this binning node and its children.
An algorithm to unfold distributions from detector to truth level.
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
Create a TH2D histogram capable to hold a covariance matrix.
virtual const Int_t * GetRowIndexArray() const
void Fatal(const char *location, const char *msgfmt,...)
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance matrix contribution from input measurement uncertainties.
virtual Double_t GetDistributionAverageBinSize(Int_t axis, Bool_t includeUnderflow, Bool_t includeOverflow) const
Get average bin size on the specified axis.
virtual Int_t ScanTau(Int_t nPoint, Double_t tauMin, Double_t tauMax, TSpline **scanResult, Int_t mode=kEScanTauRhoAvg, const char *distribution=0, const char *projectionMode=0, TGraph **lCurvePlot=0, TSpline **logTauXPlot=0, TSpline **logTauYPlot=0)
Scan a function wrt tau and determine the minimum.
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=0)
Correlated one-sigma shifts correspinding to a given systematic uncertainty.
TH1 * GetBias(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
Retrieve bias vector as a new histogram.
Base class for spline implementation containing the Draw/Paint methods.
Definition: TSpline.h:20
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4770
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
void GetBinUnderflowOverflowStatus(Int_t iBin, Int_t *uStatus, Int_t *oStatus) const
Return bit maps indicating underflow and overflow status.
#define f(i)
Definition: RSha256.hxx:104
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TH2 * GetEmatrixSysUncorr(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve covariance contribution from uncorrelated (statistical) uncertainties of the response matrix...
TH1 * GetFoldedOutput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, Bool_t addBgr=kFALSE) const
Retrieve unfolding result folded back as a new histogram.
STL namespace.
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Covariance contribution from background uncorrelated uncertainty.
TUnfoldDensity(void)
Only for use by root streamer or derived classes.
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
Short_t Abs(Short_t d)
Definition: TMathBase.h:108
TH2 * GetL(const char *histogramName, const char *histogramTitle=0, Bool_t useAxisBinning=kTRUE)
Access matrix of regularisation conditions in a new histogram.
void RegularizeDistribution(ERegMode regmode, EDensityMode densityMode, const char *distribution, const char *axisSteering)
Set up regularisation conditions.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:734
TUnfoldBinning const * GetNextNode(void) const
next sister node
Class to create third splines to interpolate knots Arbitrary conditions can be introduced for first a...
Definition: TSpline.h:191
Double_t GetRhoI(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0) const
Get global correlation coefficients, possibly cumulated over several bins.
Definition: TUnfold.cxx:3466
TH2 * GetEmatrixSysBackgroundUncorr(const char *bgrSource, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve covariance contribution from uncorrelated background uncertainties.
virtual Int_t GetDimension() const
Definition: TH1.h:277
TH1 * GetLxMinusBias(const char *histogramName, const char *histogramTitle=0)
Get regularisation conditions multiplied by result vector minus bias L(x-biasScale*biasVector).
TH2 * GetEmatrixTotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Get covariance matrix including all contributions.
const TUnfoldBinning * GetOutputBinning(const char *distributionName=0) const
Locate a binning node for the unfolded (truth level) quantities.
Double_t x[n]
Definition: legend1.C:17
static TString Format(const char *fmt,...)
Static method which formats a string using a printf style format descriptor and return a TString...
Definition: TString.cxx:2286
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
TH1 * GetBackground(const char *histogramName, const char *bgrSource=0, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, Int_t includeError=3) const
Retrieve a background source in a new histogram.
EHistMap
arrangement of axes for the response matrix (TH2 histogram)
Definition: TUnfold.h:140
EConstraint
type of extra constraint
Definition: TUnfold.h:110
Double_t Log10(Double_t x)
Definition: TMath.h:763
void RegularizeOneDistribution(const TUnfoldBinning *binning, ERegMode regmode, EDensityMode densityMode, const char *axisSteering)
Regularize the distribution of the given node.
void Info(const char *location, const char *msgfmt,...)
Int_t GetBinNeighbours(Int_t globalBin, Int_t axis, Int_t *prev, Double_t *distPrev, Int_t *next, Double_t *distNext, Bool_t isPeriodic=kFALSE) const
Get neighbour bins along the specified axis.
Int_t Finite(Double_t x)
Check if it is finite with a mask in order to be consistent in presence of fast math.
Definition: TMath.h:770
Int_t GetStartBin(void) const
first bin of this node
void Error(const char *location, const char *msgfmt,...)
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:120
Int_t GetDistributionDimension(void) const
query dimension of this node&#39;s distribution
void GetBias(TH1 *bias, const Int_t *binMap=0) const
Get bias vector including bias scale.
Definition: TUnfold.cxx:2922
void DecodeAxisSteering(const char *axisSteering, const char *options, Int_t *isOptionGiven) const
Decode axis steering.
virtual TString GetOutputBinName(Int_t iBinX) const
Get bin name of an output bin.
Definition: TUnfold.cxx:1684
ROOT::R::TRInterface & r
Definition: Object.C:4
Service class for 2-Dim histogram classes.
Definition: TH2.h:30
Class to manage histogram axis.
Definition: TAxis.h:30
TH2 * GetRhoIJtotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve correlation coefficients, including all uncertainties.
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
TH1 * GetInput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
Retrieve input distribution in a new histogram.
const TUnfoldBinning * GetInputBinning(const char *distributionName=0) const
Locate a binning node for the input (measured) quantities.
TAxis * GetYaxis()
Definition: TH1.h:316
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=0)
Correlated one-sigma shifts from shifting tau.
TH1 * GetDeltaSysTau(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve 1-sigma shift corresponding to the previously specified uncertainty on tau.
void Warning(const char *location, const char *msgfmt,...)
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
static TH2D * CreateHistogramOfMigrations(TUnfoldBinning const *xAxis, TUnfoldBinning const *yAxis, char const *histogramName, Bool_t originalXAxisBinning=kFALSE, Bool_t originalYAxisBinning=kFALSE, char const *histogramTitle=0)
Create a TH2D histogram capable to hold the bins of the two input binning schemes on the x and y axes...
Int_t GetNrows() const
Definition: TMatrixTBase.h:122
const Bool_t kFALSE
Definition: RtypesCore.h:88
TH1 * GetRhoIstatbgr(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=0)
Retrieve global correlation coefficients including input (statistical) and background uncertainties...
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
#define d(i)
Definition: RSha256.hxx:102
TH1 * GetDeltaSysSource(const char *source, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve a correlated systematic 1-sigma shift.
TString GetDistributionAxisLabel(Int_t axis) const
get name of an axis
virtual TString GetOutputBinName(Int_t iBinX) const
Get bin name of an output bin.
#define ClassImp(name)
Definition: Rtypes.h:359
double Double_t
Definition: RtypesCore.h:55
TH1 * GetOutput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
retrieve unfolding result as a new histogram
void GetProbabilityMatrix(TH2 *A, EHistMap histmap) const
Get matrix of probabilities.
Definition: TUnfold.cxx:2996
void GetFoldedOutput(TH1 *folded, const Int_t *binMap=0) const
Get unfolding result on detector level.
Definition: TUnfold.cxx:2948
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:56
void GetInput(TH1 *inputData, const Int_t *binMap=0) const
Input vector of measurements.
Definition: TUnfold.cxx:3031
Int_t GetDistributionNumberOfBins(void) const
number of bins in the distribution possibly including under/overflow
void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b, Double_t &c, Double_t &d)
Definition: TSpline.h:231
TH1 * GetDeltaSysBackgroundScale(const char *bgrSource, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Retrieve systematic 1-sigma shift corresponding to a background scale uncertainty.
void GetOutput(TH1 *output, const Int_t *binMap=0) const
Get output distribution, possibly cumulated over several bins.
Definition: TUnfold.cxx:3263
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TF1 * f1
Definition: legend1.C:11
you should not use this method at all Int_t Int_t Double_t Double_t Double_t Int_t Double_t Double_t Double_t Double_t b
Definition: TRolke.cxx:630
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:84
virtual Double_t GetScanVariable(Int_t mode, const char *distribution, const char *projectionMode)
Calculate the function for ScanTau().
#define c(i)
Definition: RSha256.hxx:101
Double_t GetDensityFactor(EDensityMode densityMode, Int_t iBin) const
Density correction factor for a given bin.
TH2 * GetEmatrixInput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Get covariance contribution from the input uncertainties (data statistical uncertainties).
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2500
TH1 * GetRhoItotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=0)
Retrieve global correlation coefficients including all uncertainty sources.
void GetL(TH2 *l) const
Get matrix of regularisation conditions.
Definition: TUnfold.cxx:3153
Definition: first.py:1
virtual Int_t GetNbinsX() const
Definition: TH1.h:291
Double_t Sqrt(Double_t x)
Definition: TMath.h:690
Int_t GetNbins() const
Definition: TAxis.h:121
void GetRhoItotal(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0)
Get global correlatiocn coefficients, summing up all contributions.
float * q
Definition: THbookFile.cxx:87
const Bool_t kTRUE
Definition: RtypesCore.h:87
Double_t Eval(Double_t x) const
Eval this spline at x.
Definition: TSpline.cxx:782
TUnfoldBinning * AddBinning(TUnfoldBinning *binning)
Add a TUnfoldBinning as the last child of this node.
const Int_t n
Definition: legend1.C:16
TUnfoldBinning const * FindNode(char const *name) const
Traverse the tree and return the first node which matches the given name.
TH2 * GetProbabilityMatrix(const char *histogramName, const char *histogramTitle=0, Bool_t useAxisBinning=kTRUE) const
Get matrix of probabilities in a new histogram.
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
TAxis * GetXaxis()
Get the behaviour adopted by the object about the statoverflows. See EStatOverflows for more informat...
Definition: TH1.h:315
virtual Int_t GetNbinsY() const
Definition: TH1.h:292
void RegularizeDistributionRecursive(const TUnfoldBinning *binning, ERegMode regmode, EDensityMode densityMode, const char *distribution, const char *axisSteering)
Recursively add regularisation conditions for this node and its children.