class TUnfold: public TObject

 TUnfold solves the inverse problem

   chi**2 = (y-Ax)# Vyy^-1 (y-Ax) + tau^2 (L(x-x0))# L(x-x0) + lambda sum_i(y_i -(Ax)_i)

 where # means that the matrix is transposed

 Monte Carlo input

   y: vector of measured quantities  (dimension ny)
   Vyy: covariance matrix for y (dimension ny x ny)
        in many cases V is diagonal and calculated from the errors of y
   A: migration matrix               (dimension ny x nx)
   x: unknown underlying distribution (dimension nx)

 Regularisation

   tau: parameter, defining the regularisation strength
   L: matrix of regularisation conditions (dimension nl x nx)
   x0: bias distribution

 Preservation of the area

   lambda: lagrangian multiplier
   y_i: one component of the vector y
   (Ax)_i: one component of the vector Ax

 and chi**2 is minimized
   (a) not constrained: minimisation is performed a function of x for fixed lambda=0
  or
   (b) constrained: minimisation is performed a function of x and lambda

 This applies to a very large number of problems, where the measured
 distribution y is a linear superposition of several Monte Carlo shapes
 and the sum of these shapes gives the output distribution x

 The constraint can be useful to reduce biases on the result x
 in cases where the vector y follows non-Gaussian probability densities
 (example: Poisson statistics at counting experiments in particle physics)

 Some random examples:

   (1) measure a cross-section as a function of, say, E_T(detector)
        and unfold it to obtain the underlying distribution E_T(generator)
   (2) measure a lifetime distribution and unfold the contributions from
        different flavours
   (3) measure the transverse mass and decay angle
        and unfold for the true mass distribution plus background

 Documentation

 Some technical documentation is available here:
   http://www.desy.de/~sschmitt

 References:

 A nice overview of the method is given in:
  The L-curve and Its Use in the Numerical Treatment of Inverse Problems
  (2000) by P. C. Hansen, in Computational Inverse Problems in
  Electrocardiology, ed. P. Johnston,
  Advances in Computational Bioengineering
  http://www.imm.dtu.dk/~pch/TR/Lcurve.ps
 The relevant equations are (1), (2) for the unfolding
 and (14) for the L-curve curvature definition

 Related literature on unfolding:
  The program package RUN and the web-page by V.Blobel
    http://www.desy.de/~blobel/unfold.html
  Talk by V. Blobel, Terascale Statistics school
    https://indico.desy.de/contributionDisplay.py?contribId=23&confId=1149
  References quoted in Blobel's talk:
   Per Chistian Hansen, Rank-Deficient and Discrete Ill-posed Problems,
        Siam (1998)
   Jari Kaipio and Erkki Somersalo, Statistical and Computational
        Inverse problems, Springer (2005)



 Implementation

 The result of the unfolding is calculated as follows:

    Lsquared = L#L            regularisation conditions squared

    epsilon_j = sum_i A_ij    vector of efficiencies

    E^-1  = ((A# Vyy^-1 A)+tau^2 Lsquared)

    x = E (A# Vyy^-1 y + tau^2 Lsquared x0 +lambda/2 * epsilon) is the result

 The derivatives
    dx_k/dy_i
    dx_k/dA_ij
    dx_k/d(tau^2)
 are calculated for further usage.

 The covariance matrix V_xx is calculated as:
    Vxx_ij = sum_kl dx_i/dy_k Vyy_kl dx_j/dy_l

 Warning:

  The algorithm is based on "standard" matrix inversion, with the
  known limitations in numerical accuracy and computing cost for
  matrices with large dimensions.

  Thus the algorithm should not used for large dimensions of x and y
    nx should not be much larger than 200
    ny should not be much larger than 1000



 Example of using TUnfold:

 imagine a 2-dimensional histogram is filled, named A
    y-axis: generated quantity (e.g. 10 bins)
    x-axis: reconstructed quantity (e.g. 20 bin)
 The data are filled in a 1-dimensional histogram, named y
 Note1: ALWAYS choose a higher number of bins on the reconstructed side
         as compared to the generated size!
 Note2: the events which are generated but not reconstructed
         have to be added to the appropriate overflow bins of A
 Note3: make sure all bins have sufficient statistics and their error is
         non-zero. By default, bins with zero error are simply skipped;
         however, this may cause problems if You try to unfold something
         which depends on these input bins.

 code fragment (with histograms A and y filled):

      TUnfold unfold(A,TUnfold::kHistMapOutputHoriz);
      Double_t tau=1.E-4;
      Double_t biasScale=0.0;
      unfold.DoUnfold(tau,y,biasScale);
      TH1D *x=unfold.GetOutput("x","myVariable");
      TH2D *rhoij=unfold.GetRhoIJ("correlation","myVariable");

 will create histograms "x" and "correlation" from A and y.
 if tau is very large, the output is biased to the generated distribution scaled by biasScale
 if tau is very small, the output will show oscillations
 and large entries in the correlation matrix


 Proper choice of tau

 One of the difficult questions is about the choice of tau. The most
 common method is the L-curve method: a two-dimensional curve is plotted
   x-axis: log10(chisquare)
   y-axis: log10(regularisation condition)
 In many cases this curve has an L-shape. The best choice of tau is in the
 kink of the L

 Within TUnfold a simple version of the L-curve analysis is included.
 It tests a given number of points in a predefined tau-range and searches
 for the maximum of the curvature in the L-curve (kink position).
 if no tau range is given, the range of teh scan is determied automatically

 Example: scan tau and produce the L-curve plot

 Code fragment: assume A and y are filled

      TUnfold unfold(A,TUnfold::kHistMapOutputHoriz);

      unfold.SetInput(y);

      Int_t nScan=30;
      Int_t iBest;
      TSpline *logTauX,*logTauY;
      TGraph *lCurve;

      iBest=unfold.ScanLcurve(nScan,0.0,0.0,&lCurve);

      std::cout<<"tau="<<unfold.GetTau()<<"\n";

      TH1D *x=unfold.GetOutput("x","myVariable");
      TH2D *rhoij=unfold.GetRhoIJ("correlation","myVariable");

 This creates
    logTauX: the L-curve's x-coordinate as a function of log(tau)
    logTauY: the L-curve's y-coordinate as a function of log(tau)
    lCurve: a graph of the L-curve
    x,rhoij: unfolding result for best choice of tau
    iBest: the coordinate/spline knot number with the best choice of tau

 Note: always check the L curve after unfolding. The algorithm is not
    perfect


 Bin averaging of the output

 Sometimes it is useful to unfold for a fine binning in x and
 calculate the final result with a smaller number of bins. The advantage
 is a reduction in the correlation coefficients if bins are averaged.
 For this type of averaging the full error matrix has to be used.
 There are methods in TUnfold to support this type of calculation
 Example:
    The vector x has dimension 49, it consists of 7x7 bins
      in two variables (Pt,Eta)
    The unfolding result is to be presented as one-dimensional projections
    in (Pt) and (Eta)
    The bins of x are mapped as: bins 1..7 the first Eta bin
                                 bins 2..14 the second Eta bin

                                 bins 1,8,15,... the first Pt bin

 code fragment:

      TUnfold unfold(A,TUnfold::kHistMapOutputHoriz);
      Double_t tau=1.E-4;
      Double_t biasScale=0.0;
      unfold.DoUnfold(tau,y,biasScale);
      Int_t binMapEta[49+2];
      Int_t binMapPt[49+2];
      // overflow and underflow bins are not used
      binMapEta[0]=-1;
      binMapEta[49+1]=-1;
      binMapPt[0]=-1;
      binMapPt[49+1]=-1;
      for(Int_t i=1;i<=49;i++) {
         // all bins (i) with the same (i-1)/7 are added
         binMapEta[i] = (i-1)/7 +1;
         // all bins (i) with the same (i-1)%7 are added
         binMapPt[i]  = (i-1)%7 +1;
      }
      TH1D *etaHist=new TH1D("eta(unfolded)",";eta",7,etamin,etamax);
      TH1D *etaCorr=new TH2D("eta(unfolded)",";eta;eta",7,etamin,etamax,7,etamin,etamax);
      TH1D *ptHist=new TH1D("pt(unfolded)",";pt",7,ptmin,ptmax);
      TH1D *ptCorr=new TH2D("pt(unfolded)",";pt;pt",7,ptmin,ptmax,7,ptmin,ptmax);
      unfold.GetOutput(etaHist,binMapEta);
      unfold.GetRhoIJ(etaCorrt,binMapEta);
      unfold.GetOutput(ptHist,binMapPt);
      unfold.GetRhoIJ(ptCorrt,binMapPt);



 Alternative Regularisation conditions

 Regularisation is needed for most unfolding problems, in order to avoid
 large oscillations and large correlations on the output bins.
 It means that some extra conditions are applied on the output bins

 Within TUnfold these conditions are posed on the difference (x-x0), where
    x:  unfolding output
    x0: the bias distribution, by default calculated from
        the input matrix A. There is a method SetBias() to change the
        bias distribution.
        The 3rd argument to DoUnfold() is a scale factor applied to the bias
          bias_default[j] = sum_i A[i][j]
          x0[j] = scaleBias*bias[j]
        The scale factor can be used to
         (a) completely suppress the bias by setting it to zero
         (b) compensate differences in the normalisation between data
             and Monte Carlo

 If the regularisation is strong, i.e. large parameter tau,
 then the distribution x or its derivatives will look like the bias
 distribution. If the parameter tau is small, the distribution x is
 independent of the bias.

 Three basic types of regularisation are implemented in TUnfold

    condition            regularisation

    kRegModeNone         none
    kRegModeSize         minimize the size of (x-x0)
    kRegModeDerivative   minimize the 1st derivative of (x-x0)
    kRegModeCurvature    minimize the 2nd derivative of (x-x0)

 kRegModeSize is the regularisation scheme which usually is found in
 literature. In addition, the bias usually is not present
 (bias scale factor is zero).

 The non-standard regularisation schemes kRegModeDerivative and
 kRegModeCurvature have the nice feature that they create correlations
 between x-bins, whereas the non-regularized unfolding tends to create
 negative correlations between bins. For these regularisation schemes the
 parameter tau could be tuned such that the correlations are smallest,
 as an alternative to the L-curve method.

 If kRegModeSize is chosen or if x is a smooth function through all bins,
 the regularisation condition can be set on all bins together by giving
 the appropriate argument in the constructor (see examples above).

 If x is composed of independent groups of bins (for example,
 signal and background binning in two variables), it may be necessary to
 set regularisation conditions for the individual groups of bins.
 In this case,  give  kRegModeNone  in the constructor and specify
 the bin grouping with calls to
          RegularizeBins()   specify a 1-dimensional group of bins
          RegularizeBins2D() specify a 2-dimensional group of bins

 For ultimate flexibility, the regularisation condition can be set on each
 bin individually
  -> give  kRegModeNone  in the constructor and use
      RegularizeSize()        regularize one bin
      RegularizeDerivative()  regularize the slope given by two bins
      RegularizeCurvature()   regularize the curvature given by three bins