Logo ROOT  
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
6An algorithm to unfold distributions from detector to truth level
7
8TUnfoldDensity is used to decompose a measurement y into several sources x,
9given the measurement uncertainties, background b and a matrix of migrations A.
10The method can be applied to a large number of problems,
11where the measured distribution y is a linear superposition
12of several Monte Carlo shapes. Beyond such a simple template fit,
13TUnfoldDensity has an adjustable regularisation term and also supports an
14optional constraint on the total number of events.
15Background sources can be specified, with a normalisation constant and
16normalisation uncertainty. In addition, variants of the response
17matrix may be specified, these are taken to determine systematic
18uncertainties. Complex, multidimensional arrangements of signal and
19background bins are managed with the help of the class TUnfoldBinning.
20
21If you use this software, please consider the following citation
22
23<b>S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]</b>
24
25Detailed documentation and updates are available on
26http://www.desy.de/~sschmitt
27
28### Brief recipe to use TUnfoldSys:
29
30 - Set up binning schemes for the truth and measured
31distributions. The binning schemes may be coded in the XML language,
32for 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
49A detailed documentation of the various GetXXX() methods to control
50systematic uncertainties is given with the method TUnfoldSys.
51
52### Why to use complex binning schemes
53
54in literature on unfolding, the "standard" test case is a
55one-dimensional distribution without underflow or overflow bins.
56The migration matrix is almost diagonal.
57
58<b>This "standard" case is rarely realized for real problems.</b>
59
60Often one has to deal with multi-dimensional distributions.
61In addition, there are underflow and overflow bins
62or other background bins, possibly determined with the help of auxiliary
63measurements.
64
65In TUnfoldDensity, such complex binning schemes are handled with the help
66of the class TUnfoldBinning. For both the measurement and the truth
67there is a tree structure. The tree nodes may correspond to single
68bins (e.g. nuisance parameters) or may hold multi-dimensional distributions.
69
70For example, the "measurement" tree could have two leaves, one for
71the primary distribution and one for auxiliary measurements.
72Similarly, the "truth" tree could have two leaves, one for the
73signal and one for the background.
74Each of the leaves may then have a multi-dimensional distribution.
75
76The class TUnfoldBinning takes care to map all bins of the
77"measurement" to a one-dimensional vector y.
78Similarly, the "truth" bins are mapped to the vector x.
79
80### How to choose the regularisation settings
81
82In 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
90Each of the algorithms has its own advantages and disadvantages.
91The algorithm (1) does not work if the input data are too similar to the
92MC 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
99The algorithm (2) only works if the variable does have a real minimum
100as a function of tau. If global correlations are minimized, the situation
101is as follows:
102The matrix of migration typically introduces negative correlations.
103The area constraint introduces some positive correlation.
104Regularisation on the "size" introduces no correlation.
105Regularisation on 1st or 2nd derivatives adds positive correlations.
106
107For these reasons, "size" regularisation does not work well with
108the tau-scan: the higher tau, the smaller rho, but there is no minimum.
109As a result, large values of tau (too strong regularisation) are found.
110In contrast, the tau-scan is expected to work better with 1st or 2nd
111derivative regularisation, because at some point the negative
112correlations from migrations are approximately cancelled by the
113positive correlations from the regularisation conditions.
114
115whichever 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 "TGraph.h"
155#include <iostream>
156#include <map>
157
158//#define DEBUG
159
160#ifdef DEBUG
161using namespace std;
162#endif
163
165
167{
168 // clean up
172}
173
174////////////////////////////////////////////////////////////////////////////////
175/// Only for use by root streamer or derived classes.
176
178{
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{
225 // set up binning schemes
226 fConstOutputBins = outputBins;
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
240 new TUnfoldBinning(*genAxis,1,1);
242 }
243 // check whether binning scheme is valid
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
254 new TUnfoldBinning(*detAxis,0,0);
256 }
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();
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();
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) {
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
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
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
432 fRegularisationConditions=new TUnfoldBinning("regularisation");
433
434 TUnfoldBinning *thisRegularisationBinning=
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#ifdef DEBUG
474 Int_t nbin=0;
475#endif
476 for(Int_t bin=startBin;bin<endBin;bin++) {
477 factor[bin-startBin]=GetDensityFactor(densityMode,bin);
478#ifdef DEBUG
479 if(factor[bin-startBin] !=0.0) nbin++;
480#endif
481 }
482#ifdef DEBUG
483 cout<<"initial number of bins "<<nbin<<"\n";
484#endif
485 Int_t dimension=binning->GetDistributionDimension();
486
487 // decide whether to skip underflow/overflow bins
488#ifdef DEBUG
489 nbin=0;
490#endif
491 for(Int_t bin=startBin;bin<endBin;bin++) {
492 Int_t uStatus,oStatus;
493 binning->GetBinUnderflowOverflowStatus(bin,&uStatus,&oStatus);
494 if(uStatus & isOptionGiven[1]) factor[bin-startBin]=0.;
495 if(oStatus & isOptionGiven[3]) factor[bin-startBin]=0.;
496#ifdef DEBUG
497 if(factor[bin-startBin] !=0.0) nbin++;
498#endif
499 }
500#ifdef DEBUG
501 cout<<"after underflow/overflow bin removal "<<nbin<<"\n";
502#endif
503 if(regmode==kRegModeSize) {
504 Int_t nRegBins=0;
505 // regularize all bins of the distribution, possibly excluding
506 // underflow/overflow bins
507 for(Int_t bin=startBin;bin<endBin;bin++) {
508 if(factor[bin-startBin]==0.0) continue;
509 if(AddRegularisationCondition(bin,factor[bin-startBin])) {
510 nRegBins++;
511 }
512 }
513 if(nRegBins) {
514 thisRegularisationBinning->AddBinning("size",nRegBins);
515 }
516 } else if((regmode==kRegModeDerivative)||(regmode==kRegModeCurvature)) {
517 for(Int_t direction=0;direction<dimension;direction++) {
518 // for each direction
519 Int_t nRegBins=0;
520 Int_t directionMask=(1<<direction);
521 if(isOptionGiven[7] & directionMask) {
522#ifdef DEBUG
523 cout<<"skip direction "<<direction<<"\n";
524#endif
525 continue;
526 }
527 Double_t binDistanceNormalisation=
528 (isOptionGiven[5] & directionMask) ?
530 (direction,isOptionGiven[0] & directionMask,
531 isOptionGiven[2] & directionMask) : 1.0;
532 for(Int_t bin=startBin;bin<endBin;bin++) {
533 // check whether bin is excluded
534 if(factor[bin-startBin]==0.0) continue;
535 // for each bin, find the neighbour bins
536 Int_t iPrev,iNext;
537 Double_t distPrev,distNext;
538 Int_t error=binning->GetBinNeighbours
539 (bin,direction,&iPrev,&distPrev,&iNext,&distNext,
540 isOptionGiven[6] & directionMask);
541 if(error) {
542 Error("RegularizeOneDistribution",
543 "invalid option %s (isPeriodic) for axis %s"
544 " (has underflow or overflow)",axisSteering,
545 binning->GetDistributionAxisLabel(direction).Data());
546 }
547 if((regmode==kRegModeDerivative)&&(iNext>=0)) {
548 Double_t f0 = -factor[bin-startBin];
549 Double_t f1 = factor[iNext-startBin];
550 if(isOptionGiven[4] & directionMask) {
551 if(distNext>0.0) {
552 f0 *= binDistanceNormalisation/distNext;
553 f1 *= binDistanceNormalisation/distNext;
554 } else {
555 f0=0.;
556 f1=0.;
557 }
558 }
559 if((f0==0.0)||(f1==0.0)) continue;
560 if(AddRegularisationCondition(bin,f0,iNext,f1)) {
561 nRegBins++;
562#ifdef DEBUG
563 std::cout<<"Added Reg: bin "<<bin<<" "<<f0
564 <<" next: "<<iNext<<" "<<f1<<"\n";
565#endif
566 }
567 } else if((regmode==kRegModeCurvature)&&(iPrev>=0)&&(iNext>=0)) {
568 Double_t f0 = factor[iPrev-startBin];
569 Double_t f1 = -factor[bin-startBin];
570 Double_t f2 = factor[iNext-startBin];
571 if(isOptionGiven[4] & directionMask) {
572 if((distPrev<0.)&&(distNext>0.)) {
573 distPrev= -distPrev;
574 Double_t f=TMath::Power(binDistanceNormalisation,2.)/
575 (distPrev+distNext);
576 f0 *= f/distPrev;
577 f1 *= f*(1./distPrev+1./distNext);
578 f2 *= f/distNext;
579 } else {
580 f0=0.;
581 f1=0.;
582 f2=0.;
583 }
584 }
585 if((f0==0.0)||(f1==0.0)||(f2==0.0)) continue;
586 if(AddRegularisationCondition(iPrev,f0,bin,f1,iNext,f2)) {
587 nRegBins++;
588#ifdef DEBUG
589 std::cout<<"Added Reg: prev "<<iPrev<<" "<<f0
590 <<" bin: "<<bin<<" "<<f1
591 <<" next: "<<iNext<<" "<<f2<<"\n";
592#endif
593 }
594 }
595 }
596 if(nRegBins) {
598 if(regmode==kRegModeDerivative) {
599 name="derivative_";
600 } else if(regmode==kRegModeCurvature) {
601 name="curvature_";
602 }
603 name += binning->GetDistributionAxisLabel(direction);
604 thisRegularisationBinning->AddBinning(name,nRegBins);
605 }
606 }
607 }
608#ifdef DEBUG
609 //fLsquared->Print();
610#endif
611}
612
613///////////////////////////////////////////////////////////////////////
614/// retrieve unfolding result as a new histogram
615///
616/// \param[in] histogramName name of the histogram
617/// \param[in] histogramTitle (default=0) title of the histogram
618/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
619/// \param[in] axisSteering (default=0) detailed steering within the extracted
620/// distribution
621/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
622/// proper binning and axis labels
623///
624/// return value: pointer to a new histogram. If
625/// <b>useAxisBinning</b> is set and if the selected distribution fits
626/// into a root histogram (1,2,3 dimensions) then return a histogram
627/// with the proper binning on each axis. Otherwise, return a 1D
628/// histogram with equidistant binning. If the histogram title is
629/// zero, a title is assigned automatically.
630///
631/// The <b>axisSteering</b> is defines as follows: "axis[mode];axis[mode];..."
632/// where:
633///
634/// - axis = name of an axis or *
635/// - mode = any combination of the letters CUO0123456789
636///
637/// - C collapse axis into one bin (add up results). If
638/// any of the numbers 0-9 are given in addition, only these bins are added up.
639/// Here bins are counted from zero, whereas in root, bins are counted
640/// from 1. Obviously, this only works for up to 10 bins.
641/// - U discard underflow bin
642/// - O discard overflow bin
643///
644/// examples: imagine the binning has two axis, named x and y.
645///
646/// - "*[UO]" exclude underflow and overflow bins for all axis.
647/// So here a TH2 is returned but all undeflow and overflow bins are empty
648/// - "x[UOC123]" integrate over the variable x but only using the
649/// bins 1,2,3 and not the underflow and overflow in x.
650/// Here a TH1 is returned, the axis is labelled "y" and
651/// the underflow and overflow (in y) are filled. However only the x-bins
652/// 1,2,3 are used to determine the content.
653/// - "x[C];y[UO]" integrate over the variable x, including
654/// underflow and overflow but exclude underflow and overflow in y.
655/// Here a TH1 is returned, the axis is labelled "y". The underflow
656/// and overflow in y are empty.
657
659(const char *histogramName,const char *histogramTitle,
660 const char *distributionName,const char *axisSteering,
661 Bool_t useAxisBinning) const
662{
663 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
664 Int_t *binMap=0;
665 TH1 *r=binning->CreateHistogram
666 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
667 if(r) {
668 TUnfoldSys::GetOutput(r,binMap);
669 }
670 if(binMap) {
671 delete [] binMap;
672 }
673 return r;
674}
675
676////////////////////////////////////////////////////////////////////////////////
677/// Retrieve bias vector as a new histogram.
678///
679/// \param[in] histogramName name of the histogram
680/// \param[in] histogramTitle (default=0) title of the histogram
681/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
682/// \param[in] axisSteering (default=0) detailed steering within the extracted
683/// distribution
684/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
685/// proper binning and axis labels
686///
687/// returns a new histogram. See method GetOutput() for a detailed
688/// description of the arguments
689
691(const char *histogramName,const char *histogramTitle,
692 const char *distributionName,const char *axisSteering,
693 Bool_t useAxisBinning) const
694{
695 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
696 Int_t *binMap=0;
697 TH1 *r=binning->CreateHistogram
698 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
699 if(r) {
700 TUnfoldSys::GetBias(r,binMap);
701 }
702 if(binMap) delete [] binMap;
703 return r;
704}
705
706////////////////////////////////////////////////////////////////////////////////
707/// Retrieve unfolding result folded back as a new histogram.
708///
709/// \param[in] histogramName name of the histogram
710/// \param[in] histogramTitle (default=0) title of the histogram
711/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
712/// \param[in] axisSteering (default=0) detailed steering within the extracted
713/// distribution
714/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
715/// proper binning and axis labels
716/// \param[in] addBgr (default=false) if true, include the background
717/// contribution (useful for direct comparison to data)
718///
719/// returns a new histogram. See method GetOutput() for a detailed
720/// description of the arguments
721
723(const char *histogramName,const char *histogramTitle,
724 const char *distributionName,const char *axisSteering,
725 Bool_t useAxisBinning,Bool_t addBgr) const
726{
727 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
728 Int_t *binMap=0;
729 TH1 *r=binning->CreateHistogram
730 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
731 if(r) {
733 if(addBgr) {
735 }
736 }
737 if(binMap) delete [] binMap;
738 return r;
739}
740
741////////////////////////////////////////////////////////////////////////////////
742/// Retrieve a background source in a new histogram.
743///
744/// \param[in] histogramName name of the histogram
745/// \param[in] bgrSource the background source to retrieve
746/// \param[in] histogramTitle (default=0) title of the histogram
747/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
748/// \param[in] axisSteering (default=0) detailed steering within the extracted
749/// distribution
750/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
751/// proper binning and axis labels
752/// \param[in] includeError (default=3) type of background errors to
753/// be included (+1 uncorrelated bgr errors, +2 correlated bgr errors)
754///
755/// returns a new histogram. See method GetOutput() for a detailed
756/// description of the arguments
757
759(const char *histogramName,const char *bgrSource,const char *histogramTitle,
760 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning,
761 Int_t includeError) const
762{
763 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
764 Int_t *binMap=0;
765 TH1 *r=binning->CreateHistogram
766 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
767 if(r) {
768 TUnfoldSys::GetBackground(r,bgrSource,binMap,includeError,kTRUE);
769 }
770 if(binMap) delete [] binMap;
771 return r;
772}
773
774////////////////////////////////////////////////////////////////////////////////
775/// Retrieve input distribution in a new histogram.
776///
777/// \param[in] histogramName name of the histogram
778/// \param[in] histogramTitle (default=0) title of the histogram
779/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
780/// \param[in] axisSteering (default=0) detailed steering within the extracted
781/// distribution
782/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
783/// proper binning and axis labels
784///
785/// returns a new histogram. See method GetOutput() for a detailed
786/// description of the arguments
787
789(const char *histogramName,const char *histogramTitle,
790 const char *distributionName,const char *axisSteering,
791 Bool_t useAxisBinning) const
792{
793 TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
794 Int_t *binMap=0;
795 TH1 *r=binning->CreateHistogram
796 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
797 if(r) {
798 TUnfoldSys::GetInput(r,binMap);
799 }
800 if(binMap) delete [] binMap;
801 return r;
802}
803
804////////////////////////////////////////////////////////////////////////////////
805/// Retrieve global correlation coefficients including all uncertainty sources.
806///
807/// \param[in] histogramName name of the histogram
808/// \param[in] histogramTitle (default=0) title of the histogram
809/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
810/// \param[in] axisSteering (default=0) detailed steering within the extracted
811/// distribution
812/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
813/// proper binning and axis labels
814/// \param[out] ematInv (default=0) to return the inverse covariance matrix
815///
816/// returns a new histogram. See method GetOutput() for a detailed
817/// description of the arguments. The inverse of the covariance matrix
818/// is stored in a new histogram returned by <b>ematInv</b> if that
819/// pointer is non-zero.
820
822(const char *histogramName,const char *histogramTitle,
823 const char *distributionName,const char *axisSteering,
824 Bool_t useAxisBinning,TH2 **ematInv) {
825 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
826 Int_t *binMap=0;
827 TH1 *r=binning->CreateHistogram
828 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
829 if(r) {
830 TH2 *invEmat=0;
831 if(ematInv) {
832 if(r->GetDimension()==1) {
833 TString ematName(histogramName);
834 ematName += "_inverseEMAT";
835 Int_t *binMap2D=0;
836 invEmat=binning->CreateErrorMatrixHistogram
837 (ematName,useAxisBinning,&binMap2D,histogramTitle,
838 axisSteering);
839 if(binMap2D) delete [] binMap2D;
840 } else {
841 Error("GetRhoItotal",
842 "can not return inverse of error matrix for this binning");
843 }
844 }
845 TUnfoldSys::GetRhoItotal(r,binMap,invEmat);
846 if(invEmat) {
847 *ematInv=invEmat;
848 }
849 }
850 if(binMap) delete [] binMap;
851 return r;
852}
853
854////////////////////////////////////////////////////////////////////////////////
855/// Retrieve global correlation coefficients including input
856/// (statistical) and background uncertainties.
857///
858/// \param[in] histogramName name of the histogram
859/// \param[in] histogramTitle (default=0) title of the histogram
860/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
861/// \param[in] axisSteering (default=0) detailed steering within the extracted
862/// distribution
863/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
864/// proper binning and axis labels
865/// \param[out] ematInv (default=0) to return the inverse covariance matrix
866///
867/// returns a new histogram. See method GetOutput() for a detailed
868/// description of the arguments. The inverse of the covariance matrix
869/// is stored in a new histogram returned by <b>ematInv</b> if that
870/// pointer is non-zero.
871
873(const char *histogramName,const char *histogramTitle,
874 const char *distributionName,const char *axisSteering,
875 Bool_t useAxisBinning,TH2 **ematInv) {
876 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
877 Int_t *binMap=0;
878 TH1 *r=binning->CreateHistogram
879 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
880 if(r) {
881 TH2 *invEmat=0;
882 if(ematInv) {
883 if(r->GetDimension()==1) {
884 TString ematName(histogramName);
885 ematName += "_inverseEMAT";
886 Int_t *binMap2D=0;
887 invEmat=binning->CreateErrorMatrixHistogram
888 (ematName,useAxisBinning,&binMap2D,histogramTitle,
889 axisSteering);
890 if(binMap2D) delete [] binMap2D;
891 } else {
892 Error("GetRhoItotal",
893 "can not return inverse of error matrix for this binning");
894 }
895 }
896 TUnfoldSys::GetRhoI(r,binMap,invEmat);
897 if(invEmat) {
898 *ematInv=invEmat;
899 }
900 }
901 if(binMap) delete [] binMap;
902 return r;
903}
904
905////////////////////////////////////////////////////////////////////////////////
906/// Retrieve a correlated systematic 1-sigma shift.
907///
908/// \param[in] source identifier of the systematic uncertainty source
909/// \param[in] histogramName name of the histogram
910/// \param[in] histogramTitle (default=0) title of the histogram
911/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
912/// \param[in] axisSteering (default=0) detailed steering within the extracted
913/// distribution
914/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
915/// proper binning and axis labels
916///
917/// returns a new histogram. See method GetOutput() for a detailed
918/// description of the arguments
919
921(const char *source,const char *histogramName,
922 const char *histogramTitle,const char *distributionName,
923 const char *axisSteering,Bool_t useAxisBinning) {
924 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
925 Int_t *binMap=0;
926 TH1 *r=binning->CreateHistogram
927 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
928 if(r) {
929 if(!TUnfoldSys::GetDeltaSysSource(r,source,binMap)) {
930 delete r;
931 r=0;
932 }
933 }
934 if(binMap) delete [] binMap;
935 return r;
936}
937
938////////////////////////////////////////////////////////////////////////////////
939/// Retrieve systematic 1-sigma shift corresponding to a background
940/// scale uncertainty.
941///
942/// \param[in] bgrSource identifier of the background
943/// \param[in] histogramName name of the histogram
944/// \param[in] histogramTitle (default=0) title of the histogram
945/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
946/// \param[in] axisSteering (default=0) detailed steering within the extracted
947/// distribution
948/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
949/// proper binning and axis labels
950///
951/// returns a new histogram. See method GetOutput() for a detailed
952/// description of the arguments
953
955(const char *bgrSource,const char *histogramName,
956 const char *histogramTitle,const char *distributionName,
957 const char *axisSteering,Bool_t useAxisBinning) {
958 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
959 Int_t *binMap=0;
960 TH1 *r=binning->CreateHistogram
961 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
962 if(r) {
963 if(!TUnfoldSys::GetDeltaSysBackgroundScale(r,bgrSource,binMap)) {
964 delete r;
965 r=0;
966 }
967 }
968 if(binMap) delete [] binMap;
969 return r;
970}
971
972////////////////////////////////////////////////////////////////////////////////
973/// Retrieve 1-sigma shift corresponding to the previously specified uncertainty
974/// on tau.
975///
976/// \param[in] histogramName name of the histogram
977/// \param[in] histogramTitle (default=0) title of the histogram
978/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
979/// \param[in] axisSteering (default=0) detailed steering within the extracted
980/// distribution
981/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
982/// proper binning and axis labels
983///
984/// returns a new histogram. See method GetOutput() for a detailed
985/// description of the arguments
986
988(const char *histogramName,const char *histogramTitle,
989 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning)
990{
991 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
992 Int_t *binMap=0;
993 TH1 *r=binning->CreateHistogram
994 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
995 if(r) {
996 if(!TUnfoldSys::GetDeltaSysTau(r,binMap)) {
997 delete r;
998 r=0;
999 }
1000 }
1001 if(binMap) delete [] binMap;
1002 return r;
1003}
1004
1005////////////////////////////////////////////////////////////////////////////////
1006/// Retrieve correlation coefficients, including all uncertainties.
1007///
1008/// \param[in] histogramName name of the histogram
1009/// \param[in] histogramTitle (default=0) title of the histogram
1010/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1011/// \param[in] axisSteering (default=0) detailed steering within the extracted
1012/// distribution
1013/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1014/// proper binning and axis labels
1015///
1016/// returns a new histogram. See method GetOutput() for a detailed
1017/// description of the arguments
1018
1020(const char *histogramName,const char *histogramTitle,
1021 const char *distributionName,const char *axisSteering,
1022 Bool_t useAxisBinning)
1023{
1025 (histogramName,histogramTitle,distributionName,
1026 axisSteering,useAxisBinning);
1027 if(r) {
1028 for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
1029 Double_t e_i=r->GetBinContent(i,i);
1030 if(e_i>0.0) e_i=TMath::Sqrt(e_i);
1031 else e_i=0.0;
1032 for(Int_t j=0;j<=r->GetNbinsY()+1;j++) {
1033 if(i==j) continue;
1034 Double_t e_j=r->GetBinContent(j,j);
1035 if(e_j>0.0) e_j=TMath::Sqrt(e_j);
1036 else e_j=0.0;
1037 Double_t e_ij=r->GetBinContent(i,j);
1038 if((e_i>0.0)&&(e_j>0.0)) {
1039 r->SetBinContent(i,j,e_ij/e_i/e_j);
1040 } else {
1041 r->SetBinContent(i,j,0.0);
1042 }
1043 }
1044 }
1045 for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
1046 if(r->GetBinContent(i,i)>0.0) {
1047 r->SetBinContent(i,i,1.0);
1048 } else {
1049 r->SetBinContent(i,i,0.0);
1050 }
1051 }
1052 }
1053 return r;
1054}
1055
1056////////////////////////////////////////////////////////////////////////////////
1057/// Retrieve covariance contribution from uncorrelated (statistical)
1058/// uncertainties of the response matrix.
1059///
1060/// \param[in] histogramName name of the histogram
1061/// \param[in] histogramTitle (default=0) title of the histogram
1062/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1063/// \param[in] axisSteering (default=0) detailed steering within the extracted
1064/// distribution
1065/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1066/// proper binning and axis labels
1067///
1068/// returns a new histogram. See method GetOutput() for a detailed
1069/// description of the arguments
1070
1072(const char *histogramName,const char *histogramTitle,
1073 const char *distributionName,const char *axisSteering,
1074 Bool_t useAxisBinning)
1075{
1076 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1077 Int_t *binMap=0;
1079 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1080 if(r) {
1082 }
1083 if(binMap) delete [] binMap;
1084 return r;
1085}
1086
1087////////////////////////////////////////////////////////////////////////////////
1088/// Retrieve covariance contribution from uncorrelated background uncertainties.
1089///
1090/// \param[in] bgrSource identifier of the background
1091/// \param[in] histogramName name of the histogram
1092/// \param[in] histogramTitle (default=0) title of the histogram
1093/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1094/// \param[in] axisSteering (default=0) detailed steering within the extracted
1095/// distribution
1096/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1097/// proper binning and axis labels
1098///
1099/// returns a new histogram. See method GetOutput() for a detailed
1100/// description of the arguments
1101
1103(const char *bgrSource,const char *histogramName,
1104 const char *histogramTitle,const char *distributionName,
1105 const char *axisSteering,Bool_t useAxisBinning)
1106{
1107 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1108 Int_t *binMap=0;
1110 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1111 if(r) {
1113 }
1114 if(binMap) delete [] binMap;
1115 return r;
1116}
1117
1118////////////////////////////////////////////////////////////////////////////////
1119/// Get covariance contribution from the input uncertainties (data
1120/// statistical uncertainties).
1121///
1122/// \param[in] histogramName name of the histogram
1123/// \param[in] histogramTitle (default=0) title of the histogram
1124/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1125/// \param[in] axisSteering (default=0) detailed steering within the extracted
1126/// distribution
1127/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1128/// proper binning and axis labels
1129///
1130/// returns a new histogram. See method GetOutput() for a detailed
1131/// description of the arguments
1132
1134(const char *histogramName,const char *histogramTitle,
1135 const char *distributionName,const char *axisSteering,
1136 Bool_t useAxisBinning)
1137{
1138 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1139 Int_t *binMap=0;
1141 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1142 if(r) {
1144 }
1145 if(binMap) delete [] binMap;
1146 return r;
1147}
1148
1149////////////////////////////////////////////////////////////////////////////////
1150/// Get matrix of probabilities in a new histogram.
1151///
1152/// \param[in] histogramName name of the histogram
1153/// \param[in] histogramTitle (default=0) title of the histogram
1154/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1155/// proper binning and axis labels
1156///
1157/// returns a new histogram. if histogramTitle is null, choose a title
1158/// automatically.
1159
1161(const char *histogramName,const char *histogramTitle,
1162 Bool_t useAxisBinning) const
1163{
1165 (fConstOutputBins,fConstInputBins,histogramName,
1166 useAxisBinning,useAxisBinning,histogramTitle);
1168 return r;
1169}
1170
1171////////////////////////////////////////////////////////////////////////////////
1172/// Get covariance matrix including all contributions.
1173///
1174/// \param[in] histogramName name of the histogram
1175/// \param[in] histogramTitle (default=0) title of the histogram
1176/// \param[in] distributionName (default=0) identifier of the distribution to be extracted
1177/// \param[in] axisSteering (default=0) detailed steering within the extracted
1178/// distribution
1179/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1180/// proper binning and axis labels
1181///
1182/// returns a new histogram. See method GetOutput() for a detailed
1183/// description of the arguments
1184
1186(const char *histogramName,const char *histogramTitle,
1187 const char *distributionName,const char *axisSteering,
1188 Bool_t useAxisBinning)
1189{
1190 TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1191 Int_t *binMap=0;
1193 (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1194 if(r) {
1196 }
1197 if(binMap) delete [] binMap;
1198 return r;
1199}
1200
1201////////////////////////////////////////////////////////////////////////////////
1202/// Access matrix of regularisation conditions in a new histogram.
1203///
1204/// \param[in] histogramName name of the histogram
1205/// \param[in] histogramTitle (default=0) title of the histogram
1206/// \param[in] useAxisBinning (default=true) if set to true, try to extract a histogram with
1207/// proper binning and axis labels
1208///
1209/// returns a new histogram. if histogramTitle is null, choose a title
1210/// automatically.
1211
1213(const char *histogramName,const char *histogramTitle,Bool_t useAxisBinning)
1214{
1218 Warning("GetL",
1219 "remove invalid scheme of regularisation conditions %d %d",
1223 }
1225 fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1226 Warning("GetL","create flat regularisation conditions scheme");
1227 }
1230 useAxisBinning,useAxisBinning,histogramTitle);
1232 return r;
1233}
1234
1235////////////////////////////////////////////////////////////////////////////////
1236/// Get regularisation conditions multiplied by result vector minus bias
1237/// L(x-biasScale*biasVector).
1238///
1239/// \param[in] histogramName name of the histogram
1240/// \param[in] histogramTitle (default=0) title of the histogram
1241///
1242/// returns a new histogram.
1243/// This is a measure of the level of regularization required per
1244/// regularisation condition.
1245/// If there are (negative or positive) spikes,
1246/// these regularisation conditions dominate
1247/// over the other regularisation conditions and may introduce
1248/// the largest biases.
1249
1251(const char *histogramName,const char *histogramTitle)
1252{
1258 Warning("GetLxMinusBias",
1259 "remove invalid scheme of regularisation conditions %d %d",
1263 }
1265 fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1266 Warning("GetLxMinusBias","create flat regularisation conditions scheme");
1267 }
1269 (histogramName,kFALSE,0,histogramTitle);
1270 const Int_t *Ldx_rows=Ldx->GetRowIndexArray();
1271 const Double_t *Ldx_data=Ldx->GetMatrixArray();
1272 for(Int_t row=0;row<Ldx->GetNrows();row++) {
1273 if(Ldx_rows[row]<Ldx_rows[row+1]) {
1274 r->SetBinContent(row+1,Ldx_data[Ldx_rows[row]]);
1275 }
1276 }
1277 delete Ldx;
1278 return r;
1279}
1280
1281////////////////////////////////////////////////////////////////////////////////
1282/// Locate a binning node for the input (measured) quantities.
1283///
1284/// \param[in] distributionName (default=0) distribution to look
1285/// for. if zero, return the root node
1286///
1287/// returns: pointer to a TUnfoldBinning object or zero if not found
1288
1290(const char *distributionName) const
1291{
1292 // find binning scheme, input bins
1293 // distributionName : the distribution to locate
1294 return fConstInputBins->FindNode(distributionName);
1295}
1296
1297////////////////////////////////////////////////////////////////////////////////
1298/// Locate a binning node for the unfolded (truth level) quantities.
1299///
1300/// \param[in] distributionName (default=0) distribution to look
1301/// for. if zero, return the root node
1302///
1303/// returns: pointer to a TUnfoldBinning object or zero if not found
1304
1306(const char *distributionName) const
1307{
1308 // find binning scheme, output bins
1309 // distributionName : the distribution to locate
1310 return fConstOutputBins->FindNode(distributionName);
1311}
1312
1313////////////////////////////////////////////////////////////////////////////////
1314/// Scan a function wrt tau and determine the minimum.
1315///
1316/// \param[in] nPoint number of points to be scanned
1317/// \param[in] tauMin smallest tau value to study
1318/// \param[in] tauMax largest tau value to study
1319/// \param[out] scanResult the scanned function wrt log(tau)
1320/// \param[in] mode 1st parameter for the scan function
1321/// \param[in] distribution 2nd parameter for the scan function
1322/// \param[in] axisSteering 3rd parameter for the scan function
1323/// \param[out] lCurvePlot for monitoring, shows the L-curve
1324/// \param[out] logTauXPlot for monitoring, L-curve(X) as a function of log(tau)
1325/// \param[out] logTauYPlot for monitoring, L-curve(Y) as a function of log(tau)
1326///
1327/// Return value: the coordinate number on the curve <b>scanResult</b>
1328/// which corresponds to the minimum
1329///
1330/// The function is scanned by repeating the following steps <b>nPoint</b>
1331/// times
1332///
1333/// 1. Choose a value of tau
1334/// 2. Perform the unfolding for this choice of tau, DoUnfold(tau)
1335/// 3. Determine the scan variable GetScanVariable()
1336///
1337/// The method GetScanVariable() defines scans of correlation
1338/// coefficients, where <b>mode</b> is chosen from the enum
1339/// EScanTauMode. In addition one may set <b>distribution</b>
1340/// and/or <b>projectionMode</b> to refine the calculation of
1341/// correlations (e.g. restrict the calculation to the signal
1342/// distribution and/or exclude underflow and overflow bins).
1343/// See the documentation of GetScanVariable() for details.
1344/// Alternative scan variables may be defined by overriding the
1345/// GetScanVariable() method.
1346///
1347/// Automatic choice of scan range: if (tauMin,tauMax) do not
1348/// correspond to a valid tau range (e.g. tauMin=tauMax=0.0) then
1349/// the tau range is determined automatically. Use with care!
1350
1352(Int_t nPoint,Double_t tauMin,Double_t tauMax,TSpline **scanResult,
1353 Int_t mode,const char *distribution,const char *axisSteering,
1354 TGraph **lCurvePlot,TSpline **logTauXPlot,TSpline **logTauYPlot)
1355{
1356 typedef std::map<Double_t,Double_t> TauScan_t;
1357 typedef std::map<Double_t,std::pair<Double_t,Double_t> > LCurve_t;
1358 TauScan_t curve;
1359 LCurve_t lcurve;
1360
1361 //==========================================================
1362 // algorithm:
1363 // (1) do the unfolding for nPoint-1 points
1364 // and store the results in the map
1365 // curve
1366 // (1a) store minimum and maximum tau to curve
1367 // (1b) insert additional points, until nPoint-1 values
1368 // have been calculated
1369 //
1370 // (2) determine the best choice of tau
1371 // do the unfolding for this point
1372 // and store the result in
1373 // curve
1374 // (3) return the result in scanResult
1375
1376 //==========================================================
1377 // (1) do the unfolding for nPoint-1 points
1378 // and store the results in
1379 // curve
1380 // (1a) store minimum and maximum tau to curve
1381
1382 if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
1383 // here no range is given, has to be determined automatically
1384 // the maximum tau is determined from the chi**2 values
1385 // observed from unfolding without regularization
1386
1387 // first unfolding, without regularisation
1388 DoUnfold(0.0);
1389
1390 // if the number of degrees of freedom is too small, create an error
1391 if(GetNdf()<=0) {
1392 Error("ScanTau","too few input bins, NDF<=0 %d",GetNdf());
1393 }
1394 Double_t X0=GetLcurveX();
1395 Double_t Y0=GetLcurveY();
1396 Double_t y0=GetScanVariable(mode,distribution,axisSteering);
1397 Info("ScanTau","logtau=-Infinity y=%lf X=%lf Y=%lf",y0,X0,Y0);
1398 {
1399 // unfolding guess maximum tau and store it
1400 Double_t logTau=
1401 0.5*(TMath::Log10(GetChi2A()+3.*TMath::Sqrt(GetNdf()+1.0))
1402 -GetLcurveY());
1403 DoUnfold(TMath::Power(10.,logTau));
1405 Fatal("ScanTau","problem (missing regularisation?) X=%f Y=%f",
1407 }
1408 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1409 curve[logTau]=y;
1410 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1411 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1413 }
1414 // minimum tau is chosen such that it is less than
1415 // 1% different from the case of no regularisation
1416 // here, several points are inserted as needed
1417 while(((int)curve.size()<nPoint-1)&&
1418 ((TMath::Abs(GetLcurveX()-X0)>0.00432)||
1419 (TMath::Abs(GetLcurveY()-Y0)>0.00432))) {
1420 Double_t logTau=(*curve.begin()).first-0.5;
1421 DoUnfold(TMath::Power(10.,logTau));
1422 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1423 curve[logTau]=y;
1424 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1425 Info("ScanTay","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1427 }
1428 } else {
1429 Double_t logTauMin=TMath::Log10(tauMin);
1430 Double_t logTauMax=TMath::Log10(tauMax);
1431 if(nPoint>1) {
1432 // insert maximum tau
1433 DoUnfold(TMath::Power(10.,logTauMax));
1434 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1435 curve[logTauMax]=y;
1436 lcurve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
1437 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMax,y,
1439 }
1440 // insert minimum tau
1441 DoUnfold(TMath::Power(10.,logTauMin));
1442 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1443 curve[logTauMin]=y;
1444 lcurve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
1445 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMin,y,
1447 }
1448
1449 //==========================================================
1450 // (1b) insert additional points, until nPoint-1 values
1451 // have been calculated
1452 while((int)curve.size()<nPoint-1) {
1453 // insert additional points
1454 // points are inserted such that the largest interval in log(tau)
1455 // is divided into two smaller intervals
1456 // however, there is a penalty term for subdividing intervals
1457 // which are far away from the minimum
1458 TauScan_t::const_iterator i0,i1;
1459 i0=curve.begin();
1460 // locate minimum
1461 Double_t logTauYMin=(*i0).first;
1462 Double_t yMin=(*i0).second;
1463 for (; i0 != curve.end(); ++i0) {
1464 if((*i0).second<yMin) {
1465 yMin=(*i0).second;
1466 logTauYMin=(*i0).first;
1467 }
1468 }
1469 // insert additional points such that large log(tau) intervals
1470 // near the minimum rho are divided into two
1471 i0=curve.begin();
1472 i1=i0;
1473 Double_t distMax=0.0;
1474 Double_t logTau=0.0;
1475 for (++i1; i1 != curve.end(); ++i1) {
1476 Double_t dist;
1477 // check size of rho interval
1478 dist=TMath::Abs((*i0).first-(*i1).first)
1479 // penalty term if distance from rhoMax is large
1480 +0.25*TMath::Power(0.5*((*i0).first+(*i1).first)-logTauYMin,2.)/
1481 ((*curve.rbegin()).first-(*curve.begin()).first)/nPoint;
1482 if((dist<=0.0)||(dist>distMax)) {
1483 distMax=dist;
1484 logTau=0.5*((*i0).first+(*i1).first);
1485 }
1486 i0=i1;
1487 }
1488 DoUnfold(TMath::Power(10.,logTau));
1489 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1490 curve[logTau]=y;
1491 lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1492 Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1494 }
1495
1496 //==========================================================
1497 // (2) determine the best choice of tau
1498 // do the unfolding for this point
1499 // and store the result in
1500 // curve
1501
1502 Double_t cTmin=0.0;
1503 {
1504 Double_t *cTi=new Double_t[curve.size()];
1505 Double_t *cCi=new Double_t[curve.size()];
1506 Int_t n=0;
1507 for (TauScan_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
1508 cTi[n]=(*i).first;
1509 cCi[n]=(*i).second;
1510 n++;
1511 }
1512 // create rho Spline
1513 TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n);
1514 // find the maximum of the curvature
1515 // if the parameter iskip is non-zero, then iskip points are
1516 // ignored when looking for the largest curvature
1517 // (there are problems with the curvature determined from the first
1518 // few points of splineX,splineY in the algorithm above)
1519 Int_t iskip=0;
1520 if(n>3) iskip=1;
1521 if(n>6) iskip=2;
1522 Double_t cCmin=cCi[iskip];
1523 cTmin=cTi[iskip];
1524 for(Int_t i=iskip;i<n-1-iskip;i++) {
1525 // find minimum on this spline section
1526 // check boundary conditions for x[i+1]
1527 Double_t xMin=cTi[i+1];
1528 Double_t yMin=cCi[i+1];
1529 if(cCi[i]<yMin) {
1530 yMin=cCi[i];
1531 xMin=cTi[i];
1532 }
1533 // find minimum for x[i]<x<x[i+1]
1534 // get spline coefficients and solve equation
1535 // derivative(x)==0
1536 Double_t x,y,b,c,d;
1537 splineC->GetCoeff(i,x,y,b,c,d);
1538 // coefficients of quadratic equation
1539 Double_t m_p_half=-c/(3.*d);
1540 Double_t q=b/(3.*d);
1541 Double_t discr=m_p_half*m_p_half-q;
1542 if(discr>=0.0) {
1543 // solution found
1544 discr=TMath::Sqrt(discr);
1545 Double_t xx;
1546 if(m_p_half>0.0) {
1547 xx = m_p_half + discr;
1548 } else {
1549 xx = m_p_half - discr;
1550 }
1551 Double_t dx=cTi[i+1]-x;
1552 // check first solution
1553 if((xx>0.0)&&(xx<dx)) {
1554 y=splineC->Eval(x+xx);
1555 if(y<yMin) {
1556 yMin=y;
1557 xMin=x+xx;
1558 }
1559 }
1560 // second solution
1561 if(xx !=0.0) {
1562 xx= q/xx;
1563 } else {
1564 xx=0.0;
1565 }
1566 // check second solution
1567 if((xx>0.0)&&(xx<dx)) {
1568 y=splineC->Eval(x+xx);
1569 if(y<yMin) {
1570 yMin=y;
1571 xMin=x+xx;
1572 }
1573 }
1574 }
1575 // check whether this local minimum is a global minimum
1576 if(yMin<cCmin) {
1577 cCmin=yMin;
1578 cTmin=xMin;
1579 }
1580 }
1581 delete splineC;
1582 delete[] cTi;
1583 delete[] cCi;
1584 }
1585 Double_t logTauFin=cTmin;
1586 DoUnfold(TMath::Power(10.,logTauFin));
1587 {
1588 Double_t y=GetScanVariable(mode,distribution,axisSteering);
1589 curve[logTauFin]=y;
1590 lcurve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
1591 Info("ScanTau","Result logtau=%lf y=%lf X=%lf Y=%lf",logTauFin,y,
1593 }
1594 //==========================================================
1595 // (3) return the result in
1596 // scanResult lCurve logTauX logTauY
1597
1598 Int_t bestChoice=-1;
1599 if(curve.size()>0) {
1600 Double_t *y=new Double_t[curve.size()];
1601 Double_t *logT=new Double_t[curve.size()];
1602 int n=0;
1603 for (TauScan_t::const_iterator i = curve.begin(); i != curve.end(); ++i) {
1604 if(logTauFin==(*i).first) {
1605 bestChoice=n;
1606 }
1607 y[n]=(*i).second;
1608 logT[n]=(*i).first;
1609 n++;
1610 }
1611 if(scanResult) {
1612 TString name;
1613 name = TString::Format("scan(%d,",mode);
1614 if(distribution) name+= distribution;
1615 name += ",";
1616 if(axisSteering) name += axisSteering;
1617 name +=")";
1618 (*scanResult)=new TSpline3(name+"%log(tau)",logT,y,n);
1619 }
1620 delete[] y;
1621 delete[] logT;
1622 }
1623 if(lcurve.size()) {
1624 Double_t *logT=new Double_t[lcurve.size()];
1625 Double_t *x=new Double_t[lcurve.size()];
1626 Double_t *y=new Double_t[lcurve.size()];
1627 Int_t n=0;
1628 for (LCurve_t::const_iterator i = lcurve.begin(); i != lcurve.end(); ++i) {
1629 logT[n]=(*i).first;
1630 x[n]=(*i).second.first;
1631 y[n]=(*i).second.second;
1632 //cout<<logT[n]<<" "<< x[n]<<" "<<y[n]<<"\n";
1633 n++;
1634 }
1635 if(lCurvePlot) {
1636 *lCurvePlot=new TGraph(n,x,y);
1637 (*lCurvePlot)->SetTitle("L curve");
1638 }
1639 if(logTauXPlot)
1640 *logTauXPlot=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
1641 if(logTauYPlot)
1642 *logTauYPlot=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
1643 delete [] y;
1644 delete [] x;
1645 delete [] logT;
1646 }
1647 return bestChoice;
1648}
1649
1650////////////////////////////////////////////////////////////////////////////////
1651/// Calculate the function for ScanTau().
1652///
1653/// \param[in] mode the variable to be calculated
1654/// \param[in] distribution distribution for which the variable
1655/// is to be calculated
1656/// \param[in] axisSteering detailed steering for selecting bins on
1657/// the axes of the distribution (see method GetRhoItotal())
1658///
1659/// return value: the scan result for the given choice of tau (for
1660/// which the unfolding was performed prior to calling this method)
1661///
1662/// In ScanTau() the unfolding is repeated for various choices of tau.
1663/// For each tau, after unfolding, GetScanVariable() is called to
1664/// determine the scan result for this choice of tau.
1665///
1666/// the following modes are implemented
1667///
1668/// - kEScanTauRhoAvg : average (stat+bgr) global correlation
1669/// - kEScanTauRhoSquaredAvg : average (stat+bgr) global correlation squared
1670/// - kEScanTauRhoMax : maximum (stat+bgr) global correlation
1671/// - kEScanTauRhoAvgSys : average (stat+bgr+sys) global correlation
1672/// - kEScanTauRhoAvgSquaredSys : average (stat+bgr+sys) global correlation squared
1673/// - kEScanTauRhoMaxSys : maximum (stat+bgr+sys) global correlation
1674
1676(Int_t mode,const char *distribution,const char *axisSteering)
1677{
1678 Double_t r=0.0;
1679 TString name("GetScanVariable(");
1680 name += TString::Format("%d,",mode);
1681 if(distribution) name += distribution;
1682 name += ",";
1683 if(axisSteering) name += axisSteering;
1684 name += ")";
1685 TH1 *rhoi=0;
1688 rhoi=GetRhoIstatbgr(name,0,distribution,axisSteering,kFALSE);
1691 rhoi=GetRhoItotal(name,0,distribution,axisSteering,kFALSE);
1692 }
1693 if(rhoi) {
1694 Double_t sum=0.0;
1695 Double_t rhoMax=0.0;
1696 Int_t n=0;
1697 for(Int_t i=0;i<=rhoi->GetNbinsX()+1;i++) {
1698 Double_t c=rhoi->GetBinContent(i);
1699 if(c>=0.) {
1700 if(c>rhoMax) rhoMax=c;
1701 sum += c;
1702 n ++;
1703 }
1704 }
1706 r=sum/n;
1707 } else if((mode==kEScanTauRhoSquareAvg)||
1709 r=sum/n;
1710 } else {
1711 r=rhoMax;
1712 }
1713 // cout<<r<<" "<<GetRhoAvg()<<" "<<GetRhoMax()<<"\n";
1714 delete rhoi;
1715 } else {
1716 Fatal("GetScanVariable","mode %d not implemented",mode);
1717 }
1718 return r;
1719}
#define d(i)
Definition: RSha256.hxx:102
#define f(i)
Definition: RSha256.hxx:104
#define c(i)
Definition: RSha256.hxx:101
const Bool_t kFALSE
Definition: RtypesCore.h:101
const Bool_t kTRUE
Definition: RtypesCore.h:100
#define ClassImp(name)
Definition: Rtypes.h:375
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t b
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t Float_t Float_t Int_t Int_t UInt_t UInt_t Rectangle_t Int_t Int_t Window_t child
Option_t Option_t TPoint TPoint const char mode
char name[80]
Definition: TGX11.cxx:110
float * q
Definition: THbookFile.cxx:89
Class to manage histogram axis.
Definition: TAxis.h:30
Int_t GetNbins() const
Definition: TAxis.h:121
A TGraph is an object made of two arrays X and Y with npoints each.
Definition: TGraph.h:41
TH1 is the base class of all histogram classes in ROOT.
Definition: TH1.h:58
TAxis * GetXaxis()
Definition: TH1.h:322
virtual Int_t GetNbinsX() const
Definition: TH1.h:295
TAxis * GetYaxis()
Definition: TH1.h:323
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:5025
Service class for 2-D histogram classes.
Definition: TH2.h:30
Int_t GetNrows() const
Definition: TMatrixTBase.h:123
const Int_t * GetRowIndexArray() const override
const Element * GetMatrixArray() const override
const char * GetName() const override
Returns name of object.
Definition: TNamed.h:47
virtual void Warning(const char *method, const char *msgfmt,...) const
Issue warning message.
Definition: TObject.cxx:955
virtual void Error(const char *method, const char *msgfmt,...) const
Issue error message.
Definition: TObject.cxx:969
virtual void Fatal(const char *method, const char *msgfmt,...) const
Issue fatal error message.
Definition: TObject.cxx:997
virtual void Info(const char *method, const char *msgfmt,...) const
Issue info message.
Definition: TObject.cxx:943
Class to create third splines to interpolate knots Arbitrary conditions can be introduced for first a...
Definition: TSpline.h:201
void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b, Double_t &c, Double_t &d) const
Definition: TSpline.h:240
Double_t Eval(Double_t x) const override
Eval this spline at x.
Definition: TSpline.cxx:786
Base class for spline implementation containing the Draw/Paint methods.
Definition: TSpline.h:31
Basic string class.
Definition: TString.h:136
const char * Data() const
Definition: TString.h:369
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:2357
Binning schemes for use with the unfolding algorithm TUnfoldDensity.
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
Create a THxx histogram capable to hold the bins of this binning node and its children.
Int_t GetDistributionDimension(void) const
query dimension of this node's distribution
TString GetBinName(Int_t iBin) const
Get the name of a bin.
Int_t GetTH1xNumberOfBins(Bool_t originalAxisBinning=kTRUE, const char *axisSteering=nullptr) const
Return the number of histogram bins required when storing this binning in a one-dimensional histogram...
Double_t GetBinSize(Int_t iBin) const
Get N-dimensional bin size.
Int_t GetDistributionNumberOfBins(void) const
number of bins in the distribution possibly including under/overflow
virtual Double_t GetBinFactor(Int_t iBin) const
Return scaling factor for the given global bin number.
Int_t GetEndBin(void) const
last+1 bin of this node (includes children)
TUnfoldBinning const * GetParentNode(void) const
mother node
static TH2D * CreateHistogramOfMigrations(TUnfoldBinning const *xAxis, TUnfoldBinning const *yAxis, char const *histogramName, Bool_t originalXAxisBinning=kFALSE, Bool_t originalYAxisBinning=kFALSE, char const *histogramTitle=nullptr)
Create a TH2D histogram capable to hold the bins of the two input binning schemes on the x and y axes...
void GetBinUnderflowOverflowStatus(Int_t iBin, Int_t *uStatus, Int_t *oStatus) const
Return bit maps indicating underflow and overflow status.
virtual Double_t GetDistributionAverageBinSize(Int_t axis, Bool_t includeUnderflow, Bool_t includeOverflow) const
Get average bin size on the specified axis.
void DecodeAxisSteering(const char *axisSteering, const char *options, Int_t *isOptionGiven) const
Decode axis steering.
TString GetDistributionAxisLabel(Int_t axis) const
get name of an axis
TUnfoldBinning * AddBinning(TUnfoldBinning *binning)
Add a TUnfoldBinning as the last child of this node.
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.
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=nullptr, const char *histogramTitle=nullptr, const char *axisSteering=nullptr) const
Create a TH2D histogram capable to hold a covariance matrix.
Int_t GetStartBin(void) const
first bin of this node
TUnfoldBinning const * FindNode(char const *name) const
Traverse the tree and return the first node which matches the given name.
TUnfoldBinning const * GetChildNode(void) const
first daughter node
An algorithm to unfold distributions from detector to truth level.
void RegularizeOneDistribution(const TUnfoldBinning *binning, ERegMode regmode, EDensityMode densityMode, const char *axisSteering)
Regularize the distribution of the given node.
@ kEScanTauRhoAvg
average global correlation coefficient (from TUnfold::GetRhoI())
@ kEScanTauRhoMax
maximum global correlation coefficient (from TUnfold::GetRhoI())
@ kEScanTauRhoSquareAvgSys
average global correlation coefficient squared (from TUnfoldSys::GetRhoItotal())
@ kEScanTauRhoMaxSys
maximum global correlation coefficient (from TUnfoldSys::GetRhoItotal())
@ kEScanTauRhoSquareAvg
average global correlation coefficient squared (from TUnfold::GetRhoI())
@ kEScanTauRhoAvgSys
average global correlation coefficient (from TUnfoldSys::GetRhoItotal())
TH1 * GetRhoItotal(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=nullptr)
Retrieve global correlation coefficients including all uncertainty sources.
TH2 * GetEmatrixSysUncorr(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve covariance contribution from uncorrelated (statistical) uncertainties of the response matrix...
Double_t GetDensityFactor(EDensityMode densityMode, Int_t iBin) const
Density correction factor for a given bin.
TH2 * GetRhoIJtotal(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve correlation coefficients, including all uncertainties.
TString GetOutputBinName(Int_t iBinX) const override
Get bin name of an output bin.
const TUnfoldBinning * fConstOutputBins
binning scheme for the output (truth level)
const TUnfoldBinning * GetOutputBinning(const char *distributionName=nullptr) const
Locate a binning node for the unfolded (truth level) quantities.
TH2 * GetEmatrixInput(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Get covariance contribution from the input uncertainties (data statistical uncertainties).
TUnfoldBinning * fRegularisationConditions
binning scheme for the regularisation conditions
TH2 * GetL(const char *histogramName, const char *histogramTitle=nullptr, Bool_t useAxisBinning=kTRUE)
Access matrix of regularisation conditions in a new histogram.
TUnfoldBinning * fOwnedOutputBins
pointer to output binning scheme if owned by this class
TH2 * GetEmatrixTotal(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Get covariance matrix including all contributions.
TUnfoldBinning * fOwnedInputBins
pointer to input binning scheme if owned by this class
TH1 * GetDeltaSysBackgroundScale(const char *bgrSource, const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve systematic 1-sigma shift corresponding to a background scale uncertainty.
void RegularizeDistribution(ERegMode regmode, EDensityMode densityMode, const char *distribution, const char *axisSteering)
Set up regularisation conditions.
TH2 * GetProbabilityMatrix(const char *histogramName, const char *histogramTitle=nullptr, Bool_t useAxisBinning=kTRUE) const
Get matrix of probabilities in a new histogram.
TH1 * GetOutput(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE) const
retrieve unfolding result as a new histogram
TH1 * GetFoldedOutput(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE, Bool_t addBgr=kFALSE) const
Retrieve unfolding result folded back as a new histogram.
const TUnfoldBinning * GetInputBinning(const char *distributionName=nullptr) const
Locate a binning node for the input (measured) quantities.
TH1 * GetInput(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE) const
Retrieve input distribution in a new histogram.
TH1 * GetRhoIstatbgr(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=nullptr)
Retrieve global correlation coefficients including input (statistical) and background uncertainties.
TH1 * GetDeltaSysSource(const char *source, const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve a correlated systematic 1-sigma shift.
TH1 * GetDeltaSysTau(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve 1-sigma shift corresponding to the previously specified uncertainty on tau.
~TUnfoldDensity(void) override
TH2 * GetEmatrixSysBackgroundUncorr(const char *bgrSource, const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE)
Retrieve covariance contribution from uncorrelated background uncertainties.
const TUnfoldBinning * fConstInputBins
binning scheme for the input (detector level)
TH1 * GetLxMinusBias(const char *histogramName, const char *histogramTitle=nullptr)
Get regularisation conditions multiplied by result vector minus bias L(x-biasScale*biasVector).
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.
TUnfoldDensity(void)
Only for use by root streamer or derived classes.
TH1 * GetBias(const char *histogramName, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE) const
Retrieve bias vector as a new histogram.
virtual Int_t ScanTau(Int_t nPoint, Double_t tauMin, Double_t tauMax, TSpline **scanResult, Int_t mode=kEScanTauRhoAvg, const char *distribution=nullptr, const char *projectionMode=nullptr, TGraph **lCurvePlot=nullptr, TSpline **logTauXPlot=nullptr, TSpline **logTauYPlot=nullptr)
Scan a function wrt tau and determine the minimum.
virtual Double_t GetScanVariable(Int_t mode, const char *distribution, const char *projectionMode)
Calculate the function for ScanTau().
EDensityMode
choice of regularisation scale factors to cinstruct the matrix L
@ kDensityModeUser
scale factors from user function in TUnfoldBinning
@ kDensityModeBinWidthAndUser
scale factors from multidimensional bin width and user function
@ kDensityModeBinWidth
scale factors from multidimensional bin width
TH1 * GetBackground(const char *histogramName, const char *bgrSource=nullptr, const char *histogramTitle=nullptr, const char *distributionName=nullptr, const char *projectionMode=nullptr, Bool_t useAxisBinning=kTRUE, Int_t includeError=3) const
Retrieve a background source in a new histogram.
An algorithm to unfold distributions from detector to truth level, with background subtraction and pr...
Definition: TUnfoldSys.h:55
void GetEmatrixTotal(TH2 *ematrix, const Int_t *binMap=nullptr)
Get total error matrix, summing up all contributions.
void GetRhoItotal(TH1 *rhoi, const Int_t *binMap=nullptr, TH2 *invEmat=nullptr)
Get global correlatiocn coefficients, summing up all contributions.
void GetBackground(TH1 *bgr, const char *bgrSource=nullptr, const Int_t *binMap=nullptr, Int_t includeError=3, Bool_t clearHist=kTRUE) const
Get background into a histogram.
Definition: TUnfoldSys.cxx:530
Bool_t GetDeltaSysBackgroundScale(TH1 *delta, const char *source, const Int_t *binMap=nullptr)
Correlated one-sigma shifts from background normalisation uncertainty.
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=nullptr)
Correlated one-sigma shifts from shifting tau.
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=nullptr)
Correlated one-sigma shifts correspinding to a given systematic uncertainty.
Definition: TUnfoldSys.cxx:999
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
Covariance matrix contribution from input measurement uncertainties.
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
Covariance contribution from background uncorrelated uncertainty.
void GetEmatrixSysUncorr(TH2 *ematrix, const Int_t *binMap=nullptr, Bool_t clearEmat=kTRUE)
Covariance contribution from uncorrelated uncertainties of the response matrix.
Definition: TUnfoldSys.cxx:731
TArrayI fHistToX
mapping of histogram bins to matrix indices
Definition: TUnfold.h:166
virtual Double_t GetLcurveY(void) const
Get value on y-axis of L-curve determined in recent unfolding.
Definition: TUnfold.cxx:3227
TMatrixDSparse * MultiplyMSparseM(const TMatrixDSparse *a, const TMatrixD *b) const
Multiply sparse matrix and a non-sparse matrix.
Definition: TUnfold.cxx:774
virtual Double_t DoUnfold(void)
Core unfolding algorithm.
Definition: TUnfold.cxx:291
TMatrixD * fX0
bias vector x0
Definition: TUnfold.h:160
void GetBias(TH1 *bias, const Int_t *binMap=nullptr) const
Get bias vector including bias scale.
Definition: TUnfold.cxx:2923
void GetProbabilityMatrix(TH2 *A, EHistMap histmap) const
Get matrix of probabilities.
Definition: TUnfold.cxx:2997
virtual TString GetOutputBinName(Int_t iBinX) const
Get bin name of an output bin.
Definition: TUnfold.cxx:1685
Double_t fBiasScale
scale factor for the bias
Definition: TUnfold.h:158
void GetFoldedOutput(TH1 *folded, const Int_t *binMap=nullptr) const
Get unfolding result on detector level.
Definition: TUnfold.cxx:2949
Bool_t AddRegularisationCondition(Int_t i0, Double_t f0, Int_t i1=-1, Double_t f1=0., Int_t i2=-1, Double_t f2=0.)
Add a row of regularisation conditions to the matrix L.
Definition: TUnfold.cxx:1937
void GetL(TH2 *l) const
Get matrix of regularisation conditions.
Definition: TUnfold.cxx:3154
const TMatrixD * GetX(void) const
vector of the unfolding result
Definition: TUnfold.h:238
EConstraint
type of extra constraint
Definition: TUnfold.h:109
virtual Double_t GetLcurveX(void) const
Get value on x-axis of L-curve determined in recent unfolding.
Definition: TUnfold.cxx:3217
Double_t GetRhoI(TH1 *rhoi, const Int_t *binMap=nullptr, TH2 *invEmat=nullptr) const
Get global correlation coefficients, possibly cumulated over several bins.
Definition: TUnfold.cxx:3467
ERegMode
choice of regularisation scheme
Definition: TUnfold.h:119
@ kRegModeNone
no regularisation, or defined later by RegularizeXXX() methods
Definition: TUnfold.h:122
@ kRegModeDerivative
regularize the 1st derivative of the output distribution
Definition: TUnfold.h:128
@ kRegModeSize
regularise the amplitude of the output distribution
Definition: TUnfold.h:125
@ kRegModeCurvature
regularize the 2nd derivative of the output distribution
Definition: TUnfold.h:131
void GetInput(TH1 *inputData, const Int_t *binMap=nullptr) const
Input vector of measurements.
Definition: TUnfold.cxx:3032
void GetOutput(TH1 *output, const Int_t *binMap=nullptr) const
Get output distribution, possibly cumulated over several bins.
Definition: TUnfold.cxx:3264
EHistMap
arrangement of axes for the response matrix (TH2 histogram)
Definition: TUnfold.h:139
@ kHistMapOutputHoriz
truth level on x-axis of the response matrix
Definition: TUnfold.h:142
Double_t GetChi2A(void) const
get χ2A contribution determined in recent unfolding
Definition: TUnfold.h:313
TMatrixDSparse * fL
regularisation conditions L
Definition: TUnfold.h:152
Int_t GetNdf(void) const
get number of degrees of freedom determined in recent unfolding
Definition: TUnfold.h:323
Double_t y[n]
Definition: legend1.C:17
Double_t x[n]
Definition: legend1.C:17
const Int_t n
Definition: legend1.C:16
TF1 * f1
Definition: legend1.C:11
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
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:768
Double_t Sqrt(Double_t x)
Returns the square root of x.
Definition: TMath.h:660
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Definition: TMath.h:719
Double_t Log10(Double_t x)
Returns the common (base-10) logarithm of x.
Definition: TMath.h:760
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
Definition: TMathBase.h:123
Definition: first.py:1
static uint64_t sum(uint64_t i)
Definition: Factory.cxx:2345