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