ROOT  6.06/09
Reference Guide
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.1, add scan type RhoSquare, small bug fixes with useAxisBinning
5 //
6 // History:
7 // Version 17.0, support for density regularisation, complex binning schemes, tau scan
8 
9 /** \class TUnfoldBinning
10  \ingroup Hist
11  TUnfold is used to decompose a measurement y into several sources x
12  given the measurement uncertainties and a matrix of migrations A
13 
14  More details are described with the documentation of TUnfold.
15 
16  For most applications, it is best to use TUnfoldDensity
17  instead of using TUnfoldSys or TUnfold
18 
19  If you use this software, please consider the following citation
20  S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]
21 
22  More documentation and updates are available on
23  http://www.desy.de/~sschmitt
24 
25  As compared to TUnfold, TUndolfDensity adds the following functionality
26  * background subtraction (see documentation of TUnfoldSys)
27  * error propagation (see documentation of TUnfoldSys)
28  * regularisation schemes respecting the bin widths
29  * support for complex, multidimensional input distributions
30 
31  Complex binning schemes are imposed on the measurements y and
32  on the result vector x with the help of the class TUnfoldBinning
33  The components of x or y are part of multi-dimensional distributions.
34  The bin widths along the relevant directions in these distributions
35  are used to calculate bin densities (number of events divided by bin width)
36  or to calculate derivatives taking into account the proper distance of
37  adjacent bin centers
38 
39  ## Complex binning schemes
40  in literature on unfolding, the "standard" test case is a
41  one-dimensional distribution without underflow or overflow bins.
42  The migration matrix is almost diagonal.
43 
44  This "standard" case is rarely realized for real problems.
45 
46  Often one has to deal with multi-dimensional input distributions.
47  In addition, there are underflow and overflow bins
48  or other background bins, possibly determined with the help of auxillary
49  measurements
50 
51  In TUnfoldDensity, such complex binning schemes are handled with the help
52  of the class TUnfoldBinning. For each vector there is a tree
53  structure. The tree nodes hold multi-dimensiopnal distributions
54 
55  For example, the "measurement" tree could have two leaves, one for
56  the primary distribution and one for auxillary measurements
57 
58  Similarly, the "truth" tree could have two leaves, one for the
59  signal and one for the background.
60 
61  each of the leaves may then have a multi-dimensional distribution.
62 
63  The class TUnfoldBinning takes care to map all bins of the
64  "measurement" to the one-dimensional vector y.
65  Similarly, the "truth" bins are mapped to the vector x.
66 
67  ## Choice of the regularisation
68  In TUnfoldDensity, two methods are implemented to determine tau**2
69  1. ScanLcurve() locate the tau where the L-curve plot has a "kink"
70  this function is implemented in the TUnfold class
71  2. ScanTau() finds the solution such that some variable
72  (e.g. global correlation coefficient) is minimized
73  this function is implemented in the TUnfoldDensity class,
74  such that the variable could be made depend on the binning scheme
75 
76  Each of the algorithms has its own advantages and disadvantages
77 
78  The algorithm (1) does not work if the input data are too similar to the
79  MC prediction, that is unfolding with tau=0 gives a least-square sum
80  of zero. Typical no-go cases of the L-curve scan are:
81  - (a) the number of measurements is too small (e.g. ny=nx)
82  - (b) the input data have no statistical fluctuations
83  [identical MC events are used to fill the matrix of migrations
84  and the vector y]
85 
86  The algorithm (2) only works if the variable does have a real minimum
87  as a function of tau.
88  If global correlations are minimized, the situation is as follows:
89  The matrix of migration typically introduces negative correlations.
90  The area constraint introduces some positive correlation.
91  Regularisation on the "size" introduces no correlation.
92  Regularisation on 1st or 2nd derivatives adds positive correlations.
93  For this reason, "size" regularisation does not work well with
94  the tau-scan: the higher tau, the smaller rho, but there is no minimum.
95  In contrast, the tau-scan is expected to work well with 1st or 2nd
96  derivative regularisation, because at some point the negative
97  correlations from migrations are approximately cancelled by the
98  positive correlations from the regularisation conditions.
99 
100  whichever algorithm is used, the output has to be checked:
101  1. The L-curve should have approximate L-shape
102  and the final choice of tau should not be at the very edge of the
103  scanned region
104  2. The scan result should have a well-defined minimum and the
105  final choice of tau should sit right in the minimum
106 */
107 
108 /*
109  This file is part of TUnfold.
110 
111  TUnfold is free software: you can redistribute it and/or modify
112  it under the terms of the GNU General Public License as published by
113  the Free Software Foundation, either version 3 of the License, or
114  (at your option) any later version.
115 
116  TUnfold is distributed in the hope that it will be useful,
117  but WITHOUT ANY WARRANTY; without even the implied warranty of
118  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
119  GNU General Public License for more details.
120 
121  You should have received a copy of the GNU General Public License
122  along with TUnfold. If not, see <http://www.gnu.org/licenses/>.
123 */
124 
125 #include "TUnfoldDensity.h"
126 #include <TMath.h>
127 #include <TVectorD.h>
128 #include <TObjString.h>
129 #include <iostream>
130 #include <map>
131 
132 // #define DEBUG
133 
135 
137 {
138  // empty constructor, for derived classes
139  fConstOutputBins=0;
140  fConstInputBins=0;
141  fOwnedOutputBins=0;
142  fOwnedInputBins=0;
143  fRegularisationConditions=0;
144 }
145 
147 {
148  // clean up
150  if(fOwnedInputBins) delete fOwnedInputBins;
152 }
153 
155 (const TH2 *hist_A, EHistMap histmap,ERegMode regmode,EConstraint constraint,
156  EDensityMode densityMode,const TUnfoldBinning *outputBins,
157  const TUnfoldBinning *inputBins,const char *regularisationDistribution,
158  const char *regularisationAxisSteering) :
159  TUnfoldSys(hist_A,histmap,kRegModeNone,constraint)
160 {
161  // set up unfolding matrix and regularisation scheme
162  // hist_A: matrix that describes the migrations
163  // histmap: mapping of the histogram axes to the unfolding output
164  // regmode: global regularisation mode
165  // constraint: type of constraint to use
166  // regularisationSteering: detailed steering for the regularisation
167  // see method RegularizeDistribution()
168  // outputBins: binning scheme of the output bins
169  // inputBins: binning scheme of the input bins
170 
171  fRegularisationConditions=0;
172  // set up binning schemes
173  fConstOutputBins = outputBins;
174  fOwnedOutputBins = 0;
175  TAxis const *genAxis,*detAxis;
176  if(histmap==kHistMapOutputHoriz) {
177  genAxis=hist_A->GetXaxis();
178  detAxis=hist_A->GetYaxis();
179  } else {
180  genAxis=hist_A->GetYaxis();
181  detAxis=hist_A->GetXaxis();
182  }
183  if(!fConstOutputBins) {
184  // underflow and overflow are included in the binning scheme
185  // They may be used on generator level
186  fOwnedOutputBins =
187  new TUnfoldBinning(*genAxis,1,1);
188  fConstOutputBins = fOwnedOutputBins;
189  }
190  // check whether binning scheme is valid
191  if(fConstOutputBins->GetParentNode()) {
192  Error("TUnfoldDensity",
193  "Invalid output binning scheme (node is not the root node)");
194  }
195  fConstInputBins = inputBins;
196  fOwnedInputBins = 0;
197  if(!fConstInputBins) {
198  // underflow and overflow are not included in the binning scheme
199  // They are still used to count events which have not been reconstructed
200  fOwnedInputBins =
201  new TUnfoldBinning(*detAxis,0,0);
202  fConstInputBins = fOwnedInputBins;
203  }
204  if(fConstInputBins->GetParentNode()) {
205  Error("TUnfoldDensity",
206  "Invalid input binning scheme (node is not the root node)");
207  }
208  // check whether binning scheme matches with the histogram
209  // in terms of total number of bins
210  Int_t nOut=genAxis->GetNbins();
211  Int_t nOutMapped=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins());
212  if(nOutMapped!= nOut) {
213  Error("TUnfoldDensity",
214  "Output binning incompatible number of bins %d!=%d",
215  nOutMapped, nOut);
216  }
217  // check whether binning scheme matches with the histogram
218  Int_t nInput=detAxis->GetNbins();
219  Int_t nInputMapped=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins());
220  if(nInputMapped!= nInput) {
221  Error("TUnfoldDensity",
222  "Input binning incompatible number of bins %d!=%d ",
223  nInputMapped, nInput);
224  }
225 
226  // report detailed list of excluded bins
227  for (Int_t ix = 0; ix <= nOut+1; ix++) {
228  if(fHistToX[ix]<0) {
229  Info("TUnfold","*NOT* unfolding bin %s",GetOutputBinName(ix).Data());
230  }
231  }
232 
233  // set up the regularisation here
234  if(regmode !=kRegModeNone) {
235  RegularizeDistribution
236  (regmode,densityMode,regularisationDistribution,
237  regularisationAxisSteering);
238  }
239 }
240 
242  if(!fConstOutputBins) return TUnfold::GetOutputBinName(iBinX);
243  else return fConstOutputBins->GetBinName(iBinX);
244 }
245 
247 (EDensityMode densityMode,Int_t iBin) const
248 {
249  // density correction factor for a given bin
250  // distributionName : name of the distribution within the output binning
251  // densityFlags : type of factor to calculate
252  // iBin : bin number
253  Double_t factor=1.0;
254  if((densityMode == kDensityModeBinWidth)||
255  (densityMode == kDensityModeBinWidthAndUser)) {
256  Double_t binSize=fConstOutputBins->GetBinSize(iBin);
257  if(binSize>0.0) factor /= binSize;
258  else factor=0.0;
259  }
260  if((densityMode == kDensityModeUser)||
261  (densityMode == kDensityModeBinWidthAndUser)) {
262  factor *= fConstOutputBins->GetBinFactor(iBin);
263  }
264  return factor;
265 }
266 
268 (ERegMode regmode,EDensityMode densityMode,const char *distribution,
269  const char *axisSteering)
270 {
271  // regularize distribution(s) using the given settings
272  // regmode: basic regularisation mode (see class TUnfold)
273  // densityMode: how to apply bin density corrections
274  // (normalisation to bin width or user factor)
275  // distribution: name of the distribiution where this regularisation
276  // is applied to (if zero, apply to all)
277  // axisSteering: regularisation steering specific to the axes
278  // The steering is defined as follows
279  // "steering1;steering2;...steeringN"
280  // each "steeringX" is defined as
281  // axisName:[options]
282  // axisName: the name of an axis where "options" applies
283  // the special name * matches all axes
284  // options: one of several character as follows
285  // u : exclude underflow bin from derivatives along this axis
286  // o : exclude overflow bin from derivatives along this axis
287  // U : exclude underflow bin
288  // O : exclude overflow bin
289  // b : use bin width for derivative calculation
290  // B : same as 'b' but in addition normalize to average bin width
291  //
292  // example: "*[UOB]" uses bin widths for derivatives and
293  // underflow/overflow bins are not regularized
294 
295  RegularizeDistributionRecursive(GetOutputBinning(),regmode,densityMode,
296  distribution,axisSteering);
297 }
298 
300 (const TUnfoldBinning *binning,ERegMode regmode,
301  EDensityMode densityMode,const char *distribution,const char *axisSteering) {
302  // recursively regularize distribution(s) using the given settings
303  // binning: distributions for this node an its children are considered
304  // regmode: basic regularisation mode (see class TUnfold)
305  // densityMode: how to apply bin density corrections
306  // (normalisation to bin withd or user factor)
307  // distribution: name of the distribiution where this regularisation
308  // is applied to (if zero, apply to all)
309  // axisSteering: regularisation steering specific to the axes
310  // (see method RegularizeDistribution())
311  if((!distribution)|| !TString(distribution).CompareTo(binning->GetName())) {
312  RegularizeOneDistribution(binning,regmode,densityMode,axisSteering);
313  }
314  for(const TUnfoldBinning *child=binning->GetChildNode();child;
315  child=child->GetNextNode()) {
316  RegularizeDistributionRecursive(child,regmode,densityMode,distribution,
317  axisSteering);
318  }
319 }
320 
322 (const TUnfoldBinning *binning,ERegMode regmode,
323  EDensityMode densityMode,const char *axisSteering)
324 {
325  // regularize the distribution in this node
326  // binning: the distributions to regularize
327  // regmode: basic regularisation mode (see class TUnfold)
328  // densityMode: how to apply bin density corrections
329  // (normalisation to bin withd or user factor)
330  // axisSteering: regularisation steering specific to the axes
331  // (see method RegularizeDistribution())
332  if(!fRegularisationConditions)
333  fRegularisationConditions=new TUnfoldBinning("regularisation");
334 
335  TUnfoldBinning *thisRegularisationBinning=
336  fRegularisationConditions->AddBinning(binning->GetName());
337 
338  // decode steering
339  Int_t isOptionGiven[6] = {0};
340  binning->DecodeAxisSteering(axisSteering,"uUoObB",isOptionGiven);
341  // U implies u
342  isOptionGiven[0] |= isOptionGiven[1];
343  // O implies o
344  isOptionGiven[2] |= isOptionGiven[3];
345  // B implies b
346  isOptionGiven[4] |= isOptionGiven[5];
347 #ifdef DEBUG
348  cout<<" "<<isOptionGiven[0]
349  <<" "<<isOptionGiven[1]
350  <<" "<<isOptionGiven[2]
351  <<" "<<isOptionGiven[3]
352  <<" "<<isOptionGiven[4]
353  <<" "<<isOptionGiven[5]
354  <<"\n";
355 #endif
356  Info("RegularizeOneDistribution","regularizing %s regMode=%d"
357  " densityMode=%d axisSteering=%s",
358  binning->GetName(),(Int_t) regmode,(Int_t)densityMode,
359  axisSteering ? axisSteering : "");
360  Int_t startBin=binning->GetStartBin();
361  Int_t endBin=startBin+ binning->GetDistributionNumberOfBins();
362  std::vector<Double_t> factor(endBin-startBin);
363  Int_t nbin=0;
364  for(Int_t bin=startBin;bin<endBin;bin++) {
365  factor[bin-startBin]=GetDensityFactor(densityMode,bin);
366  if(factor[bin-startBin] !=0.0) nbin++;
367  }
368 #ifdef DEBUG
369  cout<<"initial number of bins "<<nbin<<"\n";
370 #endif
371  Int_t dimension=binning->GetDistributionDimension();
372 
373  // decide whether to skip underflow/overflow bins
374  nbin=0;
375  for(Int_t bin=startBin;bin<endBin;bin++) {
376  Int_t uStatus,oStatus;
377  binning->GetBinUnderflowOverflowStatus(bin,&uStatus,&oStatus);
378  if(uStatus & isOptionGiven[1]) factor[bin-startBin]=0.;
379  if(oStatus & isOptionGiven[3]) factor[bin-startBin]=0.;
380  if(factor[bin-startBin] !=0.0) nbin++;
381  }
382 #ifdef DEBUG
383  cout<<"after underflow/overflow bin removal "<<nbin<<"\n";
384 #endif
385  if(regmode==kRegModeSize) {
386  Int_t nRegBins=0;
387  // regularize all bins of the distribution, possibly excluding
388  // underflow/overflow bins
389  for(Int_t bin=startBin;bin<endBin;bin++) {
390  if(factor[bin-startBin]==0.0) continue;
391  if(AddRegularisationCondition(bin,factor[bin-startBin])) {
392  nRegBins++;
393  }
394  }
395  if(nRegBins) {
396  thisRegularisationBinning->AddBinning("size",nRegBins);
397  }
398  } else if((regmode==kRegModeDerivative)||(regmode==kRegModeCurvature)) {
399  for(Int_t direction=0;direction<dimension;direction++) {
400  // for each direction
401  Int_t nRegBins=0;
402  Int_t directionMask=(1<<direction);
403  Double_t binDistanceNormalisation=
404  (isOptionGiven[5] & directionMask) ?
406  (direction,isOptionGiven[0] & directionMask,
407  isOptionGiven[2] & directionMask) : 1.0;
408  for(Int_t bin=startBin;bin<endBin;bin++) {
409  // check whether bin is excluded
410  if(factor[bin-startBin]==0.0) continue;
411  // for each bin, find the neighbour bins
412  Int_t iPrev,iNext;
413  Double_t distPrev,distNext;
414  binning->GetBinNeighbours
415  (bin,direction,&iPrev,&distPrev,&iNext,&distNext);
416  if((regmode==kRegModeDerivative)&&(iNext>=0)) {
417  Double_t f0 = -factor[bin-startBin];
418  Double_t f1 = factor[iNext-startBin];
419  if(isOptionGiven[4] & directionMask) {
420  if(distNext>0.0) {
421  f0 *= binDistanceNormalisation/distNext;
422  f1 *= binDistanceNormalisation/distNext;
423  } else {
424  f0=0.;
425  f1=0.;
426  }
427  }
428  if((f0==0.0)||(f1==0.0)) continue;
429  if(AddRegularisationCondition(bin,f0,iNext,f1)) {
430  nRegBins++;
431 #ifdef DEBUG
432  std::cout<<"Added Reg: bin "<<bin<<" "<<f0
433  <<" next: "<<iNext<<" "<<f1<<"\n";
434 #endif
435  }
436  } else if((regmode==kRegModeCurvature)&&(iPrev>=0)&&(iNext>=0)) {
437  Double_t f0 = factor[iPrev-startBin];
438  Double_t f1 = -factor[bin-startBin];
439  Double_t f2 = factor[iNext-startBin];
440  if(isOptionGiven[4] & directionMask) {
441  if((distPrev<0.)&&(distNext>0.)) {
442  distPrev= -distPrev;
443  Double_t f=TMath::Power(binDistanceNormalisation,2.)/
444  (distPrev+distNext);
445  f0 *= f/distPrev;
446  f1 *= f*(1./distPrev+1./distNext);
447  f2 *= f/distNext;
448  } else {
449  f0=0.;
450  f1=0.;
451  f2=0.;
452  }
453  }
454  if((f0==0.0)||(f1==0.0)||(f2==0.0)) continue;
455  if(AddRegularisationCondition(iPrev,f0,bin,f1,iNext,f2)) {
456  nRegBins++;
457 #ifdef DEBUG
458  std::cout<<"Added Reg: prev "<<iPrev<<" "<<f0
459  <<" bin: "<<bin<<" "<<f1
460  <<" next: "<<iNext<<" "<<f2<<"\n";
461 #endif
462  }
463  }
464  }
465  if(nRegBins) {
466  TString name;
467  if(regmode==kRegModeDerivative) {
468  name="derivative_";
469  } else if(regmode==kRegModeCurvature) {
470  name="curvature_";
471  }
472  name += binning->GetDistributionAxisLabel(direction);
473  thisRegularisationBinning->AddBinning(name,nRegBins);
474  }
475  }
476  }
477 #ifdef DEBUG
478  //fLsquared->Print();
479 #endif
480 }
481 
483 (const char *histogramName,const char *histogramTitle,
484  const char *distributionName,const char *axisSteering,
485  Bool_t useAxisBinning) const
486 {
487  // retreive unfolding result as histogram
488  // histogramName: name of the histogram
489  // histogramTitle: title of the histogram (could be zero)
490  // distributionName: for complex binning schemes specify the name
491  // of the requested distribution within the TUnfoldBinning
492  // object
493  // axisSteering:
494  // "pattern1;pattern2;...;patternN"
495  // patternI = axis[mode]
496  // axis = name or *
497  // mode = C|U|O
498  // C: collapse axis into one bin
499  // U: discarde underflow bin
500  // O: discarde overflow bin
501  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
502  Int_t *binMap=0;
503  TH1 *r=binning->CreateHistogram
504  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
505  if(r) {
506  TUnfoldSys::GetOutput(r,binMap);
507  }
508  if(binMap) {
509  delete [] binMap;
510  }
511  return r;
512 }
513 
515 (const char *histogramName,const char *histogramTitle,
516  const char *distributionName,const char *axisSteering,
517  Bool_t useAxisBinning) const
518 {
519  // retreive unfolding bias as histogram
520  // histogramName: name of the histogram
521  // histogramTitle: title of the histogram (could be zero)
522  // distributionName: for complex binning schemes specify the name
523  // of the requested distribution within the TUnfoldBinning
524  // object
525  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
526  Int_t *binMap=0;
527  TH1 *r=binning->CreateHistogram
528  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
529  if(r) {
530  TUnfoldSys::GetBias(r,binMap);
531  }
532  if(binMap) delete [] binMap;
533  return r;
534 }
535 
537 (const char *histogramName,const char *histogramTitle,
538  const char *distributionName,const char *axisSteering,
539  Bool_t useAxisBinning,Bool_t addBgr) const
540 {
541  // retreive unfolding result folded back by the matrix
542  // histogramName: name of the histogram
543  // histogramTitle: title of the histogram (could be zero)
544  // distributionName: for complex binning schemes specify the name
545  // of the requested distribution within the TUnfoldBinning
546  // object
547  // addBgr: true if the background shall be included
548  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
549  Int_t *binMap=0;
550  TH1 *r=binning->CreateHistogram
551  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
552  if(r) {
553  TUnfoldSys::GetFoldedOutput(r,binMap);
554  if(addBgr) {
555  TUnfoldSys::GetBackground(r,0,binMap,0,kFALSE);
556  }
557  }
558  if(binMap) delete [] binMap;
559  return r;
560 }
561 
563 (const char *histogramName,const char *bgrSource,const char *histogramTitle,
564  const char *distributionName,const char *axisSteering,Bool_t useAxisBinning,
565  Int_t includeError,Bool_t clearHist) const
566 {
567  // retreive a background source
568  // histogramName: name of the histogram
569  // bgrSource: name of the background source
570  // histogramTitle: title of the histogram (could be zero)
571  // distributionName: for complex binning schemes specify the name
572  // of the requested distribution within the TUnfoldBinning
573  // object
574  // include error: +1 if uncorrelated bgr errors should be included
575  // +2 if correlated bgr errors should be included
576  // clearHist: whether the histogram should be cleared
577  // if false, the background sources are added to the histogram
578  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
579  Int_t *binMap=0;
580  TH1 *r=binning->CreateHistogram
581  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
582  if(r) {
583  TUnfoldSys::GetBackground(r,bgrSource,binMap,includeError,clearHist);
584  }
585  if(binMap) delete [] binMap;
586  return r;
587 }
588 
590 (const char *histogramName,const char *histogramTitle,
591  const char *distributionName,const char *axisSteering,
592  Bool_t useAxisBinning) const
593 {
594  // retreive input distribution
595  // histogramName: name of the histogram
596  // histogramTitle: title of the histogram (could be zero)
597  // distributionName: for complex binning schemes specify the name
598  // of the requested distribution within the TUnfoldBinning
599  // object
600  TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
601  Int_t *binMap=0;
602  TH1 *r=binning->CreateHistogram
603  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
604  if(r) {
605  TUnfoldSys::GetInput(r,binMap);
606  }
607  if(binMap) delete [] binMap;
608  return r;
609 }
610 
612 (const char *histogramName,const char *histogramTitle,
613  const char *distributionName,const char *axisSteering,
614  Bool_t useAxisBinning,TH2 **ematInv) {
615  // retreive global correlation coefficients, total error
616  // and inverse of error matrix
617  // histogramName: name of the histogram
618  // histogramTitle: title of the histogram (could be zero)
619  // distributionName: for complex binning schemes specify the name
620  // of the requested distribution within the TUnfoldBinning
621  // object
622  // axisSteering:
623  // "pattern1;pattern2;...;patternN"
624  // patternI = axis[mode]
625  // axis = name or *
626  // mode = C|U|O
627  // C: collapse axis into one bin
628  // U: discarde underflow bin
629  // O: discarde overflow bin
630  // useAxisBinning: if true, try to use the axis bin widths
631  // on the output histogram
632  // ematInv: retreive inverse of error matrix
633  // if ematInv==0 the inverse is not returned
634  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
635  Int_t *binMap=0;
636  TH1 *r=binning->CreateHistogram
637  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
638  if(r) {
639  TH2 *invEmat=0;
640  if(ematInv) {
641  if(r->GetDimension()==1) {
642  TString ematName(histogramName);
643  ematName += "_inverseEMAT";
644  Int_t *binMap2D=0;
645  invEmat=binning->CreateErrorMatrixHistogram
646  (ematName,useAxisBinning,&binMap2D,histogramTitle,
647  axisSteering);
648  if(binMap2D) delete [] binMap2D;
649  } else {
650  Error("GetRhoItotal",
651  "can not return inverse of error matrix for this binning");
652  }
653  }
654  TUnfoldSys::GetRhoItotal(r,binMap,invEmat);
655  if(invEmat) {
656  *ematInv=invEmat;
657  }
658  }
659  if(binMap) delete [] binMap;
660  return r;
661 }
662 
664 (const char *histogramName,const char *histogramTitle,
665  const char *distributionName,const char *axisSteering,
666  Bool_t useAxisBinning,TH2 **ematInv) {
667  // retreive global correlation coefficients, input error
668  // and inverse of corresponding error matrix
669  // histogramName: name of the histogram
670  // histogramTitle: title of the histogram (could be zero)
671  // distributionName: for complex binning schemes specify the name
672  // of the requested distribution within the TUnfoldBinning
673  // object
674  // axisSteering:
675  // "pattern1;pattern2;...;patternN"
676  // patternI = axis[mode]
677  // axis = name or *
678  // mode = C|U|O
679  // C: collapse axis into one bin
680  // U: discarde underflow bin
681  // O: discarde overflow bin
682  // useAxisBinning: if true, try to use the axis bin widths
683  // on the output histogram
684  // ematInv: retreive inverse of error matrix
685  // if ematInv==0 the inverse is not returned
686  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
687  Int_t *binMap=0;
688  TH1 *r=binning->CreateHistogram
689  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
690  if(r) {
691  TH2 *invEmat=0;
692  if(ematInv) {
693  if(r->GetDimension()==1) {
694  TString ematName(histogramName);
695  ematName += "_inverseEMAT";
696  Int_t *binMap2D=0;
697  invEmat=binning->CreateErrorMatrixHistogram
698  (ematName,useAxisBinning,&binMap2D,histogramTitle,
699  axisSteering);
700  if(binMap2D) delete [] binMap2D;
701  } else {
702  Error("GetRhoItotal",
703  "can not return inverse of error matrix for this binning");
704  }
705  }
706  TUnfoldSys::GetRhoI(r,binMap,invEmat);
707  if(invEmat) {
708  *ematInv=invEmat;
709  }
710  }
711  if(binMap) delete [] binMap;
712  return r;
713 }
714 
716 (const char *source,const char *histogramName,
717  const char *histogramTitle,const char *distributionName,
718  const char *axisSteering,Bool_t useAxisBinning) {
719  // retreive histogram of systematic 1-sigma shifts
720  // source: name of systematic error
721  // histogramName: name of the histogram
722  // histogramTitle: title of the histogram (could be zero)
723  // distributionName: for complex binning schemes specify the name
724  // of the requested distribution within the TUnfoldBinning
725  // object
726  // axisSteering:
727  // "pattern1;pattern2;...;patternN"
728  // patternI = axis[mode]
729  // axis = name or *
730  // mode = C|U|O
731  // C: collapse axis into one bin
732  // U: discarde underflow bin
733  // O: discarde overflow bin
734  // useAxisBinning: if true, try to use the axis bin widths
735  // on the output histogram
736  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
737  Int_t *binMap=0;
738  TH1 *r=binning->CreateHistogram
739  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
740  if(r) {
741  if(!TUnfoldSys::GetDeltaSysSource(r,source,binMap)) {
742  delete r;
743  r=0;
744  }
745  }
746  if(binMap) delete [] binMap;
747  return r;
748 }
749 
751 (const char *bgrSource,const char *histogramName,
752  const char *histogramTitle,const char *distributionName,
753  const char *axisSteering,Bool_t useAxisBinning) {
754  // retreive histogram of systematic 1-sigma shifts due to a background
755  // normalisation uncertainty
756  // source: name of background source
757  // histogramName: name of the histogram
758  // histogramTitle: title of the histogram (could be zero)
759  // distributionName: for complex binning schemes specify the name
760  // of the requested distribution within the TUnfoldBinning
761  // object
762  // axisSteering:
763  // "pattern1;pattern2;...;patternN"
764  // patternI = axis[mode]
765  // axis = name or *
766  // mode = C|U|O
767  // C: collapse axis into one bin
768  // U: discarde underflow bin
769  // O: discarde overflow bin
770  // useAxisBinning: if true, try to use the axis bin widths
771  // on the output histogram
772  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
773  Int_t *binMap=0;
774  TH1 *r=binning->CreateHistogram
775  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
776  if(r) {
777  if(!TUnfoldSys::GetDeltaSysBackgroundScale(r,bgrSource,binMap)) {
778  delete r;
779  r=0;
780  }
781  }
782  if(binMap) delete [] binMap;
783  return r;
784 }
785 
787 (const char *histogramName,const char *histogramTitle,
788  const char *distributionName,const char *axisSteering,Bool_t useAxisBinning)
789 {
790  // retreive histogram of systematic 1-sigma shifts due to tau uncertainty
791  // histogramName: name of the histogram
792  // histogramTitle: title of the histogram (could be zero)
793  // distributionName: for complex binning schemes specify the name
794  // of the requested distribution within the TUnfoldBinning
795  // object
796  // axisSteering:
797  // "pattern1;pattern2;...;patternN"
798  // patternI = axis[mode]
799  // axis = name or *
800  // mode = C|U|O
801  // C: collapse axis into one bin
802  // U: discarde underflow bin
803  // O: discarde overflow bin
804  // useAxisBinning: if true, try to use the axis bin widths
805  // on the output histogram
806  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
807  Int_t *binMap=0;
808  TH1 *r=binning->CreateHistogram
809  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
810  if(r) {
811  if(!TUnfoldSys::GetDeltaSysTau(r,binMap)) {
812  delete r;
813  r=0;
814  }
815  }
816  if(binMap) delete [] binMap;
817  return r;
818 }
819 
821 (const char *histogramName,const char *histogramTitle,
822  const char *distributionName,const char *axisSteering,
823  Bool_t useAxisBinning)
824 {
825  // retreive histogram of total corelation coefficients, including systematic
826  // uncertainties
827  // histogramName: name of the histogram
828  // histogramTitle: title of the histogram (could be zero)
829  // distributionName: for complex binning schemes specify the name
830  // of the requested distribution within the TUnfoldBinning
831  // object
832  // axisSteering:
833  // "pattern1;pattern2;...;patternN"
834  // patternI = axis[mode]
835  // axis = name or *
836  // mode = C|U|O
837  // C: collapse axis into one bin
838  // U: discarde underflow bin
839  // O: discarde overflow bin
840  // useAxisBinning: if true, try to use the axis bin widths
841  // on the output histogram
842  TH2 *r=GetEmatrixTotal
843  (histogramName,histogramTitle,distributionName,
844  axisSteering,useAxisBinning);
845  if(r) {
846  for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
847  Double_t e_i=r->GetBinContent(i,i);
848  if(e_i>0.0) e_i=TMath::Sqrt(e_i);
849  else e_i=0.0;
850  for(Int_t j=0;j<=r->GetNbinsY()+1;j++) {
851  if(i==j) continue;
852  Double_t e_j=r->GetBinContent(j,j);
853  if(e_j>0.0) e_j=TMath::Sqrt(e_j);
854  else e_j=0.0;
855  Double_t e_ij=r->GetBinContent(i,j);
856  if((e_i>0.0)&&(e_j>0.0)) {
857  r->SetBinContent(i,j,e_ij/e_i/e_j);
858  } else {
859  r->SetBinContent(i,j,0.0);
860  }
861  }
862  }
863  for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
864  if(r->GetBinContent(i,i)>0.0) {
865  r->SetBinContent(i,i,1.0);
866  } else {
867  r->SetBinContent(i,i,0.0);
868  }
869  }
870  }
871  return r;
872 }
873 
875 (const char *histogramName,const char *histogramTitle,
876  const char *distributionName,const char *axisSteering,
877  Bool_t useAxisBinning)
878 {
879  // get error matrix contribution from uncorrelated errors on the matrix A
880  // histogramName: name of the histogram
881  // histogramTitle: title of the histogram (could be zero)
882  // distributionName: for complex binning schemes specify the name
883  // of the requested distribution within the TUnfoldBinning
884  // object
885  // axisSteering:
886  // "pattern1;pattern2;...;patternN"
887  // patternI = axis[mode]
888  // axis = name or *
889  // mode = C|U|O
890  // C: collapse axis into one bin
891  // U: discarde underflow bin
892  // O: discarde overflow bin
893  // useAxisBinning: if true, try to use the axis bin widths
894  // on the output histogram
895  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
896  Int_t *binMap=0;
898  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
899  if(r) {
901  }
902  if(binMap) delete [] binMap;
903  return r;
904 }
905 
906 
908 (const char *bgrSource,const char *histogramName,
909  const char *histogramTitle,const char *distributionName,
910  const char *axisSteering,Bool_t useAxisBinning)
911 {
912  // get error matrix from uncorrelated error of one background source
913  // histogramName: name of the histogram
914  // histogramTitle: title of the histogram (could be zero)
915  // distributionName: for complex binning schemes specify the name
916  // of the requested distribution within the TUnfoldBinning
917  // object
918  // axisSteering:
919  // "pattern1;pattern2;...;patternN"
920  // patternI = axis[mode]
921  // axis = name or *
922  // mode = C|U|O
923  // C: collapse axis into one bin
924  // U: discarde underflow bin
925  // O: discarde overflow bin
926  // useAxisBinning: if true, try to use the axis bin widths
927  // on the output histogram
928  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
929  Int_t *binMap=0;
931  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
932  if(r) {
934  }
935  if(binMap) delete [] binMap;
936  return r;
937 }
938 
940 (const char *histogramName,const char *histogramTitle,
941  const char *distributionName,const char *axisSteering,
942  Bool_t useAxisBinning)
943 {
944  // get error contribution from input vector
945  // histogramName: name of the histogram
946  // histogramTitle: title of the histogram (could be zero)
947  // distributionName: for complex binning schemes specify the name
948  // of the requested distribution within the TUnfoldBinning
949  // object
950  // axisSteering:
951  // "pattern1;pattern2;...;patternN"
952  // patternI = axis[mode]
953  // axis = name or *
954  // mode = C|U|O
955  // C: collapse axis into one bin
956  // U: discarde underflow bin
957  // O: discarde overflow bin
958  // useAxisBinning: if true, try to use the axis bin widths
959  // on the output histogram
960  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
961  Int_t *binMap=0;
963  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
964  if(r) {
965  TUnfoldSys::GetEmatrixInput(r,binMap);
966  }
967  if(binMap) delete [] binMap;
968  return r;
969 }
970 
972 (const char *histogramName,const char *histogramTitle,
973  Bool_t useAxisBinning) const
974 {
975  // get matrix of probabilities
976  // histogramName: name of the histogram
977  // histogramTitle: title of the histogram (could be zero)
978  // useAxisBinning: if true, try to get the histogram using
979  // the original matrix binning
981  (fConstOutputBins,fConstInputBins,histogramName,
982  useAxisBinning,useAxisBinning,histogramTitle);
983  TUnfold::GetProbabilityMatrix(r,kHistMapOutputHoriz);
984  return r;
985 }
986 
988 (const char *histogramName,const char *histogramTitle,
989  const char *distributionName,const char *axisSteering,
990  Bool_t useAxisBinning)
991 {
992  // get total error including systematic,statistical,background,tau errors
993  // histogramName: name of the histogram
994  // histogramTitle: title of the histogram (could be zero)
995  // distributionName: for complex binning schemes specify the name
996  // of the requested distribution within the TUnfoldBinning
997  // object
998  // axisSteering:
999  // "pattern1;pattern2;...;patternN"
1000  // patternI = axis[mode]
1001  // axis = name or *
1002  // mode = C|U|O
1003  // C: collapse axis into one bin
1004  // U: discarde underflow bin
1005  // O: discarde overflow bin
1006  // useAxisBinning: if true, try to use the axis bin widths
1007  // on the output histogram
1008  TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
1009  Int_t *binMap=0;
1010  TH2 *r=binning->CreateErrorMatrixHistogram
1011  (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
1012  if(r) {
1013  TUnfoldSys::GetEmatrixTotal(r,binMap);
1014  }
1015  if(binMap) delete [] binMap;
1016  return r;
1017 }
1018 
1020 (const char *histogramName,const char *histogramTitle,Bool_t useAxisBinning)
1021 {
1022  // return the matrix of regularisation conditions in a histogram
1023  // input:
1024  // histogramName: name of the histogram
1025  // histogramTitle: title of the histogram (could be zero)
1026  // useAxisBinning: if true, try to use the axis bin widths
1027  // on the x-axis of the output histogram
1028  if(fRegularisationConditions &&
1029  (fRegularisationConditions->GetEndBin()-
1030  fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
1031  Warning("GetL",
1032  "remove invalid scheme of regularisation conditions %d %d",
1033  fRegularisationConditions->GetEndBin(),fL->GetNrows());
1034  delete fRegularisationConditions;
1035  fRegularisationConditions=0;
1036  }
1037  if(!fRegularisationConditions) {
1038  fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1039  Warning("GetL","create flat regularisation conditions scheme");
1040  }
1042  (fConstOutputBins,fRegularisationConditions,histogramName,
1043  useAxisBinning,useAxisBinning,histogramTitle);
1044  TUnfold::GetL(r);
1045  return r;
1046 }
1047 
1049 (const char *histogramName,const char *histogramTitle)
1050 {
1051  // get regularisation conditions multiplied by result vector minus bias
1052  // L(x-biasScale*biasVector)
1053  // this is a measure of the level of regulartisation required per
1054  // regularisation condition
1055  // if there are (negative or positive) spikes,
1056  // these regularisation conditions dominate
1057  // over the other regularisation conditions
1058  // input
1059  // histogramName: name of the histogram
1060  // histogramTitle: title of the histogram (could be zero)
1061  TMatrixD dx(*GetX(), TMatrixD::kMinus, fBiasScale * (*fX0));
1062  TMatrixDSparse *Ldx=MultiplyMSparseM(fL,&dx);
1063  if(fRegularisationConditions &&
1064  (fRegularisationConditions->GetEndBin()-
1065  fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
1066  Warning("GetLxMinusBias",
1067  "remove invalid scheme of regularisation conditions %d %d",
1068  fRegularisationConditions->GetEndBin(),fL->GetNrows());
1069  delete fRegularisationConditions;
1070  fRegularisationConditions=0;
1071  }
1072  if(!fRegularisationConditions) {
1073  fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
1074  Warning("GetLxMinusBias","create flat regularisation conditions scheme");
1075  }
1076  TH1 *r=fRegularisationConditions->CreateHistogram
1077  (histogramName,kFALSE,0,histogramTitle);
1078  const Int_t *Ldx_rows=Ldx->GetRowIndexArray();
1079  const Double_t *Ldx_data=Ldx->GetMatrixArray();
1080  for(Int_t row=0;row<Ldx->GetNrows();row++) {
1081  if(Ldx_rows[row]<Ldx_rows[row+1]) {
1082  r->SetBinContent(row+1,Ldx_data[Ldx_rows[row]]);
1083  }
1084  }
1085  delete Ldx;
1086  return r;
1087 }
1088 
1090 (const char *distributionName) const
1091 {
1092  // find binning scheme, input bins
1093  // distributionName : the distribution to locate
1094  return fConstInputBins->FindNode(distributionName);
1095 }
1096 
1098 (const char *distributionName) const
1099 {
1100  // find binning scheme, output bins
1101  // distributionName : the distribution to locate
1102  return fConstOutputBins->FindNode(distributionName);
1103 }
1104 
1106 (Int_t nPoint,Double_t tauMin,Double_t tauMax,TSpline **scanResult,
1107  Int_t mode,const char *distribution,const char *axisSteering,
1108  TGraph **lCurvePlot,TSpline **logTauXPlot,TSpline **logTauYPlot)
1109 {
1110  // scan some variable as a function of tau and determine the minimum
1111  // input:
1112  // nPoint: number of points to be scanned on the resulting curve
1113  // tauMin: smallest tau value to study
1114  // tauMax: largest tau value to study
1115  // if (mauMin,tauMax) do not correspond to a valid tau range
1116  // (e.g. tauMin=tauMax=0.0) then the tau range is determined
1117  // automatically
1118  // mode,distribution,axisSteering: argument to GetScanVariable()
1119  // output:
1120  // scanResult: output spline of the variable as a function of tau
1121  // the following plots are produced on request (if pointers are non-zero)
1122  // lCurvePlot: for monitoring: the L-curve
1123  // logTauXPlot: for monitoring: L-curve(x) as a function of log(tau)
1124  // logTauYPlot: for monitoring: L-curve(y) as a function of log(tau)
1125  // return value: the coordinate number (0..nPoint-1) corresponding to the
1126  // final choice of tau
1127  typedef std::map<Double_t,Double_t> TauScan_t;
1128  typedef std::map<Double_t,std::pair<Double_t,Double_t> > LCurve_t;
1129  TauScan_t curve;
1130  LCurve_t lcurve;
1131 
1132  //==========================================================
1133  // algorithm:
1134  // (1) do the unfolding for nPoint-1 points
1135  // and store the results in the map
1136  // curve
1137  // (1a) store minimum and maximum tau to curve
1138  // (1b) insert additional points, until nPoint-1 values
1139  // have been calculated
1140  //
1141  // (2) determine the best choice of tau
1142  // do the unfolding for this point
1143  // and store the result in
1144  // curve
1145  // (3) return the result in scanResult
1146 
1147  //==========================================================
1148  // (1) do the unfolding for nPoint-1 points
1149  // and store the results in
1150  // curve
1151  // (1a) store minimum and maximum tau to curve
1152 
1153  if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
1154  // here no range is given, has to be determined automatically
1155  // the maximum tau is determined from the chi**2 values
1156  // observed from unfolding without regulatisation
1157 
1158  // first unfolding, without regularisation
1159  DoUnfold(0.0);
1160 
1161  // if the number of degrees of freedom is too small, create an error
1162  if(GetNdf()<=0) {
1163  Error("ScanTau","too few input bins, NDF<=0 %d",GetNdf());
1164  }
1165  Double_t X0=GetLcurveX();
1166  Double_t Y0=GetLcurveY();
1167  Double_t y0=GetScanVariable(mode,distribution,axisSteering);
1168  Info("ScanTau","logtau=-Infinity y=%lf X=%lf Y=%lf",y0,X0,Y0);
1169  {
1170  // unfolding guess maximum tau and store it
1171  Double_t logTau=
1172  0.5*(TMath::Log10(GetChi2A()+3.*TMath::Sqrt(GetNdf()+1.0))
1173  -GetLcurveY());
1174  DoUnfold(TMath::Power(10.,logTau));
1175  if((!TMath::Finite(GetLcurveX())) ||(!TMath::Finite(GetLcurveY()))) {
1176  Fatal("ScanTau","problem (missing regularisation?) X=%f Y=%f",
1177  GetLcurveX(),GetLcurveY());
1178  }
1179  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1180  curve[logTau]=y;
1181  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1182  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1183  GetLcurveX(),GetLcurveY());
1184  }
1185  // minimum tau is chosen such that it is less than
1186  // 1% different from the case of no regularisation
1187  // here, several points are inserted as needed
1188  while(((int)curve.size()<nPoint-1)&&
1189  ((TMath::Abs(GetLcurveX()-X0)>0.00432)||
1190  (TMath::Abs(GetLcurveY()-Y0)>0.00432))) {
1191  Double_t logTau=(*curve.begin()).first-0.5;
1192  DoUnfold(TMath::Power(10.,logTau));
1193  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1194  curve[logTau]=y;
1195  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1196  Info("ScanTay","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1197  GetLcurveX(),GetLcurveY());
1198  }
1199  } else {
1200  Double_t logTauMin=TMath::Log10(tauMin);
1201  Double_t logTauMax=TMath::Log10(tauMax);
1202  if(nPoint>1) {
1203  // insert maximum tau
1204  DoUnfold(TMath::Power(10.,logTauMax));
1205  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1206  curve[logTauMax]=y;
1207  lcurve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
1208  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMax,y,
1209  GetLcurveX(),GetLcurveY());
1210  }
1211  // insert minimum tau
1212  DoUnfold(TMath::Power(10.,logTauMin));
1213  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1214  curve[logTauMin]=y;
1215  lcurve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
1216  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMin,y,
1217  GetLcurveX(),GetLcurveY());
1218  }
1219 
1220  //==========================================================
1221  // (1b) insert additional points, until nPoint-1 values
1222  // have been calculated
1223  while((int)curve.size()<nPoint-1) {
1224  // insert additional points
1225  // points are inserted such that the largest interval in log(tau)
1226  // is divided into two smaller intervals
1227  // however, there is a penalty term for subdividing intervals
1228  // which are far away from the minimum
1229  TauScan_t::const_iterator i0,i1;
1230  i0=curve.begin();
1231  // locate minimum
1232  Double_t logTauYMin=(*i0).first;
1233  Double_t yMin=(*i0).second;
1234  for(;i0!=curve.end();i0++) {
1235  if((*i0).second<yMin) {
1236  yMin=(*i0).second;
1237  logTauYMin=(*i0).first;
1238  }
1239  }
1240  // insert additional points such that large log(tau) intervals
1241  // near the minimum rho are divided into two
1242  i0=curve.begin();
1243  i1=i0;
1244  Double_t distMax=0.0;
1245  Double_t logTau=0.0;
1246  for(i1++;i1!=curve.end();i1++) {
1247  Double_t dist;
1248  // check size of rho interval
1249  dist=TMath::Abs((*i0).first-(*i1).first)
1250  // penalty term if distance from rhoMax is large
1251  +0.25*TMath::Power(0.5*((*i0).first+(*i1).first)-logTauYMin,2.)/
1252  ((*curve.rbegin()).first-(*curve.begin()).first)/nPoint;
1253  if((dist<=0.0)||(dist>distMax)) {
1254  distMax=dist;
1255  logTau=0.5*((*i0).first+(*i1).first);
1256  }
1257  i0=i1;
1258  }
1259  DoUnfold(TMath::Power(10.,logTau));
1260  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1261  curve[logTau]=y;
1262  lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
1263  Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
1264  GetLcurveX(),GetLcurveY());
1265  }
1266 
1267  //==========================================================
1268  // (2) determine the best choice of tau
1269  // do the unfolding for this point
1270  // and store the result in
1271  // curve
1272 
1273  Double_t cTmin=0.0;
1274  {
1275  Double_t *cTi=new Double_t[curve.size()];
1276  Double_t *cCi=new Double_t[curve.size()];
1277  Int_t n=0;
1278  for(TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
1279  cTi[n]=(*i).first;
1280  cCi[n]=(*i).second;
1281  n++;
1282  }
1283  // create rho Spline
1284  TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n);
1285  // find the maximum of the curvature
1286  // if the parameter iskip is non-zero, then iskip points are
1287  // ignored when looking for the largest curvature
1288  // (there are problems with the curvature determined from the first
1289  // few points of splineX,splineY in the algorithm above)
1290  Int_t iskip=0;
1291  if(n>3) iskip=1;
1292  if(n>6) iskip=2;
1293  Double_t cCmin=cCi[iskip];
1294  cTmin=cTi[iskip];
1295  for(Int_t i=iskip;i<n-1-iskip;i++) {
1296  // find minimum on this spline section
1297  // check boundary conditions for x[i+1]
1298  Double_t xMin=cTi[i+1];
1299  Double_t yMin=cCi[i+1];
1300  if(cCi[i]<yMin) {
1301  yMin=cCi[i];
1302  xMin=cTi[i];
1303  }
1304  // find minimum for x[i]<x<x[i+1]
1305  // get spline coefficiencts and solve equation
1306  // derivative(x)==0
1307  Double_t x,y,b,c,d;
1308  splineC->GetCoeff(i,x,y,b,c,d);
1309  // coefficiencts of quadratic equation
1310  Double_t m_p_half=-c/(3.*d);
1311  Double_t q=b/(3.*d);
1312  Double_t discr=m_p_half*m_p_half-q;
1313  if(discr>=0.0) {
1314  // solution found
1315  discr=TMath::Sqrt(discr);
1316  Double_t xx;
1317  if(m_p_half>0.0) {
1318  xx = m_p_half + discr;
1319  } else {
1320  xx = m_p_half - discr;
1321  }
1322  Double_t dx=cTi[i+1]-x;
1323  // check first solution
1324  if((xx>0.0)&&(xx<dx)) {
1325  y=splineC->Eval(x+xx);
1326  if(y<yMin) {
1327  yMin=y;
1328  xMin=x+xx;
1329  }
1330  }
1331  // second solution
1332  if(xx !=0.0) {
1333  xx= q/xx;
1334  } else {
1335  xx=0.0;
1336  }
1337  // check second solution
1338  if((xx>0.0)&&(xx<dx)) {
1339  y=splineC->Eval(x+xx);
1340  if(y<yMin) {
1341  yMin=y;
1342  xMin=x+xx;
1343  }
1344  }
1345  }
1346  // check whether this local minimum is a global minimum
1347  if(yMin<cCmin) {
1348  cCmin=yMin;
1349  cTmin=xMin;
1350  }
1351  }
1352  delete splineC;
1353  delete[] cTi;
1354  delete[] cCi;
1355  }
1356  Double_t logTauFin=cTmin;
1357  DoUnfold(TMath::Power(10.,logTauFin));
1358  {
1359  Double_t y=GetScanVariable(mode,distribution,axisSteering);
1360  curve[logTauFin]=y;
1361  lcurve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
1362  Info("ScanTau","Result logtau=%lf y=%lf X=%lf Y=%lf",logTauFin,y,
1363  GetLcurveX(),GetLcurveY());
1364  }
1365  //==========================================================
1366  // (3) return the result in
1367  // scanResult lCurve logTauX logTauY
1368 
1369  Int_t bestChoice=-1;
1370  if(curve.size()>0) {
1371  Double_t *y=new Double_t[curve.size()];
1372  Double_t *logT=new Double_t[curve.size()];
1373  int n=0;
1374  for( TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
1375  if(logTauFin==(*i).first) {
1376  bestChoice=n;
1377  }
1378  y[n]=(*i).second;
1379  logT[n]=(*i).first;
1380  n++;
1381  }
1382  if(scanResult) {
1383  TString name;
1384  name = TString::Format("scan(%d,",mode);
1385  if(distribution) name+= distribution;
1386  name += ",";
1387  if(axisSteering) name += axisSteering;
1388  name +=")";
1389  (*scanResult)=new TSpline3(name+"%log(tau)",logT,y,n);
1390  }
1391  delete[] y;
1392  delete[] logT;
1393  }
1394  if(lcurve.size()) {
1395  Double_t *logT=new Double_t[lcurve.size()];
1396  Double_t *x=new Double_t[lcurve.size()];
1397  Double_t *y=new Double_t[lcurve.size()];
1398  Int_t n=0;
1399  for(LCurve_t::const_iterator i=lcurve.begin();i!=lcurve.end();i++) {
1400  logT[n]=(*i).first;
1401  x[n]=(*i).second.first;
1402  y[n]=(*i).second.second;
1403  //cout<<logT[n]<<" "<< x[n]<<" "<<y[n]<<"\n";
1404  n++;
1405  }
1406  if(lCurvePlot) {
1407  *lCurvePlot=new TGraph(n,x,y);
1408  (*lCurvePlot)->SetTitle("L curve");
1409  }
1410  if(logTauXPlot)
1411  *logTauXPlot=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
1412  if(logTauYPlot)
1413  *logTauYPlot=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
1414  delete [] y;
1415  delete [] x;
1416  delete [] logT;
1417  }
1418  return bestChoice;
1419 }
1420 
1422 (Int_t mode,const char *distribution,const char *axisSteering)
1423 {
1424  // calculate variable for ScanTau()
1425  // the unfolding is repeated for various choices of tau.
1426  // For each tau, after unfolding, the ScanTau() method calls
1427  // GetScanVariable() to determine the value of the variable which
1428  // is to be scanned
1429  //
1430  // the variable is expected to have a minimum near the "optimal" choice
1431  // of tau
1432  //
1433  // input:
1434  // mode : define the type of variable to be calculated
1435  // distribution : define the distribution for which the variable
1436  // is to be calculated
1437  // the following modes are implemented:
1438  // kEScanTauRhoAvg : average global correlation from input data
1439  // kEScanTauRhoSquaredAvg : average global correlation squared
1440  // from input data
1441  // kEScanTauRhoMax : maximum global correlation from input data
1442  // kEScanTauRhoAvgSys : average global correlation
1443  // including systematic uncertainties
1444  // kEScanTauRhoAvgSquaredSys : average global correlation squared
1445  // including systematic uncertainties
1446  // kEScanTauRhoMaxSys : maximum global correlation
1447  // including systematic uncertainties
1448  // distribution : name of the TUnfoldBinning node
1449  // for which to calculate the correlations
1450  // axisSteering : axis steering for calculating the correlations
1451  // the distribution
1452  // return: the value of the variable as determined from the present
1453  // unfolding
1454 
1455  Double_t r=0.0;
1456  TString name("GetScanVariable(");
1457  name += TString::Format("%d,",mode);
1458  if(distribution) name += distribution;
1459  name += ",";
1460  if(axisSteering) name += axisSteering;
1461  name += ")";
1462  TH1 *rhoi=0;
1463  if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoMax)||
1464  (mode==kEScanTauRhoSquareAvg)) {
1465  rhoi=GetRhoIstatbgr(name,0,distribution,axisSteering,kFALSE);
1466  } else if((mode==kEScanTauRhoAvgSys)||(mode==kEScanTauRhoMaxSys)||
1467  (mode==kEScanTauRhoSquareAvgSys)) {
1468  rhoi=GetRhoItotal(name,0,distribution,axisSteering,kFALSE);
1469  }
1470  if(rhoi) {
1471  Double_t sum=0.0;
1472  Double_t sumSquare=0.0;
1473  Double_t rhoMax=0.0;
1474  Int_t n=0;
1475  for(Int_t i=0;i<=rhoi->GetNbinsX()+1;i++) {
1476  Double_t c=rhoi->GetBinContent(i);
1477  if(c>=0.) {
1478  if(c>rhoMax) rhoMax=c;
1479  sum += c;
1480  sumSquare += c*c;
1481  n ++;
1482  }
1483  }
1484  if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoAvgSys)) {
1485  r=sum/n;
1486  } else if((mode==kEScanTauRhoSquareAvg)||
1487  (mode==kEScanTauRhoSquareAvgSys)) {
1488  r=sum/n;
1489  } else {
1490  r=rhoMax;
1491  }
1492  // cout<<r<<" "<<GetRhoAvg()<<" "<<GetRhoMax()<<"\n";
1493  delete rhoi;
1494  } else {
1495  Fatal("GetScanVariable","mode %d not implemented",mode);
1496  }
1497  return r;
1498 }
TString GetBinName(Int_t iBin) const
double dist(Rotation3D const &r1, Rotation3D const &r2)
Definition: 3DDistances.cxx:48
TH1 * GetFoldedOutput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, Bool_t addBgr=kFALSE) const
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH1.cxx:4629
const TUnfoldBinning * fConstOutputBins
void GetOutput(TH1 *output, const Int_t *binMap=0) const
Definition: TUnfold.cxx:3033
void Fatal(const char *location, const char *msgfmt,...)
void GetEmatrixInput(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
virtual const Element * GetMatrixArray() const
virtual Int_t ScanTau(Int_t nPoint, Double_t tauMin, Double_t tauMax, TSpline **scanResult, Int_t mode=kEScanTauRhoAvg, const char *distribution=0, const char *projectionMode=0, TGraph **lCurvePlot=0, TSpline **logTauXPlot=0, TSpline **logTauYPlot=0)
Bool_t GetDeltaSysSource(TH1 *hist_delta, const char *source, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:940
virtual Int_t GetDimension() const
Definition: TH1.h:283
Base class for spline implementation containing the Draw/Paint methods //.
Definition: TSpline.h:22
TH1 * GetBias(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
virtual TString GetOutputBinName(Int_t iBinX) const
Definition: TUnfold.cxx:1737
const TUnfoldBinning * GetOutputBinning(const char *distributionName=0) const
TUnfold is used to decompose a measurement y into several sources x given the measurement uncertainti...
Definition: TUnfoldSys.h:51
Int_t GetStartBin(void) const
Basic string class.
Definition: TString.h:137
int Int_t
Definition: RtypesCore.h:41
bool Bool_t
Definition: RtypesCore.h:59
TH1 * GetOutput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
const Bool_t kFALSE
Definition: Rtypes.h:92
TH2 * GetEmatrixSysUncorr(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
virtual Int_t GetNbinsX() const
Definition: TH1.h:296
void GetEmatrixSysBackgroundUncorr(TH2 *ematrix, const char *source, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
TH1 * GetInput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE) const
TUnfoldDensity(void)
TUnfoldBinning * fOwnedOutputBins
const TUnfoldBinning * GetInputBinning(const char *distributionName=0) const
Short_t Abs(Short_t d)
Definition: TMathBase.h:110
TH2 * GetL(const char *histogramName, const char *histogramTitle=0, Bool_t useAxisBinning=kTRUE)
void RegularizeDistribution(ERegMode regmode, EDensityMode densityMode, const char *distribution, const char *axisSteering)
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:501
TH2 * GetEmatrixSysBackgroundUncorr(const char *bgrSource, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
TUnfoldBinning const * GetChildNode(void) const
TH1 * GetLxMinusBias(const char *histogramName, const char *histogramTitle=0)
TH2 * GetEmatrixTotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
Double_t x[n]
Definition: legend1.C:17
void GetFoldedOutput(TH1 *folded, const Int_t *binMap=0) const
Definition: TUnfold.cxx:2815
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:2334
void GetEmatrixSysUncorr(TH2 *ematrix, const Int_t *binMap=0, Bool_t clearEmat=kTRUE)
Definition: TUnfoldSys.cxx:704
TH1 * CreateHistogram(const char *histogramName, Bool_t originalAxisBinning=kFALSE, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
EHistMap
Definition: TUnfold.h:188
Double_t Eval(Double_t x) const
Eval this spline at x.
Definition: TSpline.cxx:818
EConstraint
Definition: TUnfold.h:103
Double_t Log10(Double_t x)
Definition: TMath.h:529
void RegularizeOneDistribution(const TUnfoldBinning *binning, ERegMode regmode, EDensityMode densityMode, const char *axisSteering)
void Info(const char *location, const char *msgfmt,...)
std::vector< std::vector< double > > Data
Int_t Finite(Double_t x)
Definition: TMath.h:532
ClassImp(TUnfoldDensity) TUnfoldDensity
void Error(const char *location, const char *msgfmt,...)
TH2 * GetProbabilityMatrix(const char *histogramName, const char *histogramTitle=0, Bool_t useAxisBinning=kTRUE) const
ERegMode
Definition: TUnfold.h:107
Int_t GetNrows() const
Definition: TMatrixTBase.h:134
virtual Double_t GetBinContent(Int_t bin) const
Return content of bin number bin.
Definition: TH2.h:90
Double_t GetDensityFactor(EDensityMode densityMode, Int_t iBin) const
TUnfoldBinning * fRegularisationConditions
void GetBinUnderflowOverflowStatus(Int_t iBin, Int_t *uStatus, Int_t *oStatus) const
Double_t GetRhoI(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0) const
Definition: TUnfold.cxx:3216
ROOT::R::TRInterface & r
Definition: Object.C:4
virtual const Int_t * GetRowIndexArray() const
Service class for 2-Dim histogram classes.
Definition: TH2.h:36
Class to manage histogram axis.
Definition: TAxis.h:36
Int_t GetNbins() const
Definition: TAxis.h:125
TH2 * GetRhoIJtotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
TH1 * GetBackground(const char *histogramName, const char *bgrSource=0, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, Int_t includeError=3, Bool_t clearHist=kTRUE) const
void GetEmatrixTotal(TH2 *ematrix, const Int_t *binMap=0)
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content see convention for numbering bins in TH1::GetBin In case the bin number is greater th...
Definition: TH1.cxx:8543
TUnfoldBinning const * GetNextNode(void) const
return fString CompareTo(((TObjString *) obj) ->fString)
virtual const char * GetName() const
Returns name of object.
Definition: TNamed.h:51
Int_t GetDistributionNumberOfBins(void) const
TAxis * GetYaxis()
Definition: TH1.h:320
Bool_t GetDeltaSysTau(TH1 *delta, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:980
TH1 * GetDeltaSysTau(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
void Warning(const char *location, const char *msgfmt,...)
static TH2D * CreateHistogramOfMigrations(TUnfoldBinning const *xAxis, TUnfoldBinning const *yAxis, char const *histogramName, Bool_t originalXAxisBinning=kFALSE, Bool_t originalYAxisBinning=kFALSE, char const *histogramTitle=0)
TH1 * GetRhoIstatbgr(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=0)
TH2D * CreateErrorMatrixHistogram(const char *histogramName, Bool_t originalAxisBinning, Int_t **binMap=0, const char *histogramTitle=0, const char *axisSteering=0) const
This class serves as a container of analysis bins analysis bins are specified by defining the axes of...
TH1 * GetDeltaSysSource(const char *source, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
double f(double x)
double Double_t
Definition: RtypesCore.h:55
Int_t GetDistributionDimension(void) const
void GetBinNeighbours(Int_t globalBin, Int_t axis, Int_t *prev, Double_t *distPrev, Int_t *next, Double_t *distNext) const
virtual ~TUnfoldDensity(void)
Double_t y[n]
Definition: legend1.C:17
The TH1 histogram class.
Definition: TH1.h:80
#define name(a, b)
Definition: linkTestLib0.cpp:5
TString GetDistributionAxisLabel(Int_t axis) const
void GetCoeff(Int_t i, Double_t &x, Double_t &y, Double_t &b, Double_t &c, Double_t &d)
Definition: TSpline.h:233
virtual Int_t GetNbinsY() const
Definition: TH1.h:297
TUnfoldBinning const * FindNode(char const *name) const
TH1 * GetDeltaSysBackgroundScale(const char *bgrSource, const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
TUnfoldBinning * fOwnedInputBins
double f2(const double *x)
void GetL(TH2 *l) const
Definition: TUnfold.cxx:2977
A Graph is a graphics object made of two arrays X and Y with npoints each.
Definition: TGraph.h:53
TF1 * f1
Definition: legend1.C:11
virtual Double_t GetScanVariable(Int_t mode, const char *distribution, const char *projectionMode)
TH2 * GetEmatrixInput(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE)
virtual void SetBinContent(Int_t bin, Double_t content)
Set bin content.
Definition: TH2.cxx:2689
TH1 * GetRhoItotal(const char *histogramName, const char *histogramTitle=0, const char *distributionName=0, const char *projectionMode=0, Bool_t useAxisBinning=kTRUE, TH2 **ematInv=0)
Double_t Sqrt(Double_t x)
Definition: TMath.h:464
void GetRhoItotal(TH1 *rhoi, const Int_t *binMap=0, TH2 *invEmat=0)
virtual TString GetOutputBinName(Int_t iBinX) const
void GetBias(TH1 *bias, const Int_t *binMap=0) const
Definition: TUnfold.cxx:2794
void DecodeAxisSteering(const char *axisSteering, const char *options, Int_t *isOptionGiven) const
float * q
Definition: THbookFile.cxx:87
TUnfoldBinning * AddBinning(TUnfoldBinning *binning)
const Int_t n
Definition: legend1.C:16
void GetProbabilityMatrix(TH2 *A, EHistMap histmap) const
Definition: TUnfold.cxx:2863
void GetInput(TH1 *inputData, const Int_t *binMap=0) const
Definition: TUnfold.cxx:2884
Bool_t GetDeltaSysBackgroundScale(TH1 *delta, const char *source, const Int_t *binMap=0)
Definition: TUnfoldSys.cxx:958
virtual Double_t GetDistributionAverageBinSize(Int_t axis, Bool_t includeUnderflow, Bool_t includeOverflow) const
TAxis * GetXaxis()
Definition: TH1.h:319
void GetBackground(TH1 *bgr, const char *bgrSource=0, const Int_t *binMap=0, Int_t includeError=3, Bool_t clearHist=kTRUE) const
Definition: TUnfoldSys.cxx:480
void RegularizeDistributionRecursive(const TUnfoldBinning *binning, ERegMode regmode, EDensityMode densityMode, const char *distribution, const char *axisSteering)