// Author: Stefan Schmitt, Amnon Harel
// DESY and CERN, 11/08/11

//  Version 17.1, add scan type RhoSquare, small bug fixes with useAxisBinning
//
//  History:
//    Version 17.0, support for density regularisation, complex binning schemes, tau scan

//////////////////////////////////////////////////////////////////////////
//
//  TUnfoldDensity : public TUnfoldSys : public TUnfold
//
//  TUnfold is used to decompose a measurement y into several sources x
//  given the measurement uncertainties and a matrix of migrations A
//
//  More details are described with the documentation of TUnfold.
//
//  For most applications, it is best to use TUnfoldDensity
//  instead of using TUnfoldSys or TUnfold
//
//  If you use this software, please consider the following citation
//       S.Schmitt, JINST 7 (2012) T10003 [arXiv:1205.6201]
//
//  More documentation and updates are available on
//      http://www.desy.de/~sschmitt
//
//
//  As compared to TUnfold, TUndolfDensity adds the following functionality
//    * background subtraction (see documentation of TUnfoldSys)
//    * error propagation (see documentation of TUnfoldSys)
//    * regularisation schemes respecting the bin widths
//    * support for complex, multidimensional input distributions
//
//  Complex binning schemes are imposed on the measurements y and
//  on the result vector x with the help of the class TUnfoldBinning
//  The components of x or y are part of multi-dimensional distributions.
//  The bin widths along the relevant directions in these distributions
//  are used to calculate bin densities (number of events divided by bin width)
//  or to calculate derivatives taking into account the proper distance of 
//  adjacent bin centers
//
//  Complex binning schemes
//  =======================
//  in literature on unfolding, the "standard" test case is a
//  one-dimensional distribution without underflow or overflow bins.
//  The migration matrix is almost diagonal.
//
//  This "standard" case is rarely realized for real problems.
//
//  Often one has to deal with multi-dimensional input distributions.
//  In addition, there are underflow and overflow bins
//  or other background bins, possibly determined with the help of auxillary
//  measurements
//
//  In TUnfoldDensity, such complex binning schemes are handled with the help
//  of the class TUnfoldBinning. For each vector there is a tree
//  structure. The tree nodes hold multi-dimensiopnal distributions 
//
//  For example, the "measurement" tree could have two leaves, one for
//  the primary distribution and one for auxillary measurements
//
//  Similarly, the "truth" tree could have two leaves, one for the
//  signal and one for the background.
//
//  each of the leaves may then have a multi-dimensional distribution.
//
//  The class TUnfoldBinning takes care to map all bins of the
//  "measurement" to the one-dimensional vector y.
//  Similarly, the "truth" bins are mapped to the vector x.
//
//  Choice of the regularisation
//  ============================
//  In TUnfoldDensity, two methods are implemented to determine tau**2
//    (1)  ScanLcurve()  locate the tau where the L-curve plot has a "kink"
//      this function is implemented in the TUnfold class
//    (2)  ScanTau() finds the solution such that some variable
//           (e.g. global correlation coefficient) is minimized
//      this function is implemented in the TUnfoldDensity class,
//      such that the variable could be made depend on the binning scheme
//
//  Each of the algorithms has its own advantages and disadvantages
//
//  The algorithm (1) does not work if the input data are too similar to the
//  MC prediction, that is unfolding with tau=0 gives a least-square sum
//  of zero. Typical no-go cases of the L-curve scan are:
//    (a) the number of measurements is too small (e.g. ny=nx)
//    (b) the input data have no statistical fluctuations
//         [identical MC events are used to fill the matrix of migrations
//          and the vector y]
//
//  The algorithm (2) only works if the variable does have a real minimum
//  as a function of tau.
//  If global correlations are minimized, the situation is as follows:
//  The matrix of migration typically introduces negative correlations.
//   The area constraint introduces some positive correlation.
//   Regularisation on the "size" introduces no correlation.
//   Regularisation on 1st or 2nd derivatives adds positive correlations.
//   For this reason, "size" regularisation does not work well with
//   the tau-scan: the higher tau, the smaller rho, but there is no minimum.
//   In contrast, the tau-scan is expected to work well with 1st or 2nd
//   derivative regularisation, because at some point the negative
//   correlations from migrations are approximately cancelled by the
//   positive correlations from the regularisation conditions.
//
//  whichever algorithm is used, the output has to be checked:
//  (1) The L-curve should have approximate L-shape
//       and the final choice of tau should not be at the very edge of the
//       scanned region
//  (2) The scan result should have a well-defined minimum and the
//       final choice of tau should sit right in the minimum
//
////////////////////////////////////////////////////////////////////////////

/*
  This file is part of TUnfold.

  TUnfold is free software: you can redistribute it and/or modify
  it under the terms of the GNU General Public License as published by
  the Free Software Foundation, either version 3 of the License, or
  (at your option) any later version.

  TUnfold is distributed in the hope that it will be useful,
  but WITHOUT ANY WARRANTY; without even the implied warranty of
  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  GNU General Public License for more details.

  You should have received a copy of the GNU General Public License
  along with TUnfold.  If not, see <http://www.gnu.org/licenses/>.
*/

#include "TUnfoldDensity.h"
#include <TMath.h>
#include <TVectorD.h>
#include <TObjString.h>
#include <iostream>
#include <map>

// #define DEBUG

ClassImp(TUnfoldDensity)

TUnfoldDensity::TUnfoldDensity(void)
{
   // empty constructor, for derived classes
   fConstOutputBins=0;
   fConstInputBins=0;
   fOwnedOutputBins=0;
   fOwnedInputBins=0;
   fRegularisationConditions=0;
}

TUnfoldDensity::~TUnfoldDensity(void)
{
   // clean up
   if(fOwnedOutputBins) delete fOwnedOutputBins;
   if(fOwnedInputBins) delete fOwnedInputBins;
   if(fRegularisationConditions) delete fRegularisationConditions;
}

TUnfoldDensity::TUnfoldDensity
(const TH2 *hist_A, EHistMap histmap,ERegMode regmode,EConstraint constraint,
 EDensityMode densityMode,const TUnfoldBinning *outputBins,
 const TUnfoldBinning *inputBins,const char *regularisationDistribution,
 const char *regularisationAxisSteering) :
   TUnfoldSys(hist_A,histmap,kRegModeNone,constraint) 
{
   // set up unfolding matrix and regularisation scheme
   //    hist_A:  matrix that describes the migrations
   //    histmap: mapping of the histogram axes to the unfolding output 
   //    regmode: global regularisation mode
   //    constraint: type of constraint to use
   //    regularisationSteering: detailed steering for the regularisation
   //                  see method RegularizeDistribution()
   //    outputBins: binning scheme of the output bins
   //    inputBins: binning scheme of the input bins

   fRegularisationConditions=0;
   // set up binning schemes
   fConstOutputBins = outputBins;
   fOwnedOutputBins = 0;
   TAxis const *genAxis,*detAxis;
   if(histmap==kHistMapOutputHoriz) {
      genAxis=hist_A->GetXaxis();
      detAxis=hist_A->GetYaxis();
   } else {
      genAxis=hist_A->GetYaxis();
      detAxis=hist_A->GetXaxis();
   }
   if(!fConstOutputBins) {
      // underflow and overflow are included in the binning scheme
      // They may be used on generator level
      fOwnedOutputBins =
         new TUnfoldBinning(*genAxis,1,1);
      fConstOutputBins = fOwnedOutputBins;
   }
   // check whether binning scheme is valid
   if(fConstOutputBins->GetParentNode()) {
      Error("TUnfoldDensity",
            "Invalid output binning scheme (node is not the root node)");
   }
   fConstInputBins = inputBins;
   fOwnedInputBins = 0;
   if(!fConstInputBins) {
      // underflow and overflow are not included in the binning scheme
      // They are still used to count events which have not been reconstructed
      fOwnedInputBins =
         new TUnfoldBinning(*detAxis,0,0);
      fConstInputBins = fOwnedInputBins;
   }
   if(fConstInputBins->GetParentNode()) {
      Error("TUnfoldDensity",
            "Invalid input binning scheme (node is not the root node)");
   }
   // check whether binning scheme matches with the histogram
   // in terms of total number of bins
   Int_t nOut=genAxis->GetNbins();
   Int_t nOutMapped=TMath::Abs(fConstOutputBins->GetTH1xNumberOfBins());
   if(nOutMapped!= nOut) {
      Error("TUnfoldDensity",
            "Output binning incompatible number of bins %d!=%d",
            nOutMapped, nOut);
   }
   // check whether binning scheme matches with the histogram
   Int_t nInput=detAxis->GetNbins();
   Int_t nInputMapped=TMath::Abs(fConstInputBins->GetTH1xNumberOfBins());
   if(nInputMapped!= nInput) {
      Error("TUnfoldDensity",
            "Input binning incompatible number of bins %d!=%d ",
            nInputMapped, nInput);
   }

   // report detailed list of excluded bins
   for (Int_t ix = 0; ix <= nOut+1; ix++) {
      if(fHistToX[ix]<0) {
         Info("TUnfold","*NOT* unfolding bin %s",GetOutputBinName(ix).Data());
      }
   }

   // set up the regularisation here
   if(regmode !=kRegModeNone) {
      RegularizeDistribution
      (regmode,densityMode,regularisationDistribution,
       regularisationAxisSteering);
   }
}

TString TUnfoldDensity::GetOutputBinName(Int_t iBinX) const {
   if(!fConstOutputBins) return TUnfold::GetOutputBinName(iBinX);
   else return fConstOutputBins->GetBinName(iBinX);
}

Double_t TUnfoldDensity::GetDensityFactor
(EDensityMode densityMode,Int_t iBin) const
{
   // density correction factor for a given bin
   //    distributionName : name of the distribution within the output binning
   //    densityFlags : type of factor to calculate
   //    iBin : bin number
   Double_t factor=1.0;
   if((densityMode == kDensityModeBinWidth)||
      (densityMode == kDensityModeBinWidthAndUser)) {
      Double_t binSize=fConstOutputBins->GetBinSize(iBin);
      if(binSize>0.0) factor /= binSize;
      else factor=0.0;
   }
   if((densityMode == kDensityModeUser)||
      (densityMode == kDensityModeBinWidthAndUser)) {
      factor *= fConstOutputBins->GetBinFactor(iBin);
   }
   return factor;
}

void TUnfoldDensity::RegularizeDistribution
(ERegMode regmode,EDensityMode densityMode,const char *distribution,
 const char *axisSteering)
{
   // regularize distribution(s) using the given settings
   //     regmode: basic regularisation mode (see class TUnfold)
   //     densityMode: how to apply bin density corrections
   //              (normalisation to bin width or user factor)
   //     distribution: name of the distribiution where this regularisation
   //             is applied to (if zero, apply to all)
   //     axisSteering: regularisation steering specific to the axes
   //          The steering is defined as follows
   //             "steering1;steering2;...steeringN"
   //          each "steeringX" is defined as
   //             axisName:[options]
   //          axisName: the name of an axis where "options" applies
   //                    the special name * matches all axes
   //          options: one of several character as follows
   //             u : exclude underflow bin from derivatives along this axis
   //             o : exclude overflow bin from derivatives along this axis
   //             U : exclude underflow bin
   //             O : exclude overflow bin
   //             b : use bin width for derivative calculation
   //             B : same as 'b' but in addition normalize to average bin width
   //
   //          example:  "*[UOB]" uses bin widths for derivatives and
   //                             underflow/overflow bins are not regularized
  
   RegularizeDistributionRecursive(GetOutputBinning(),regmode,densityMode,
                                   distribution,axisSteering);
}

void TUnfoldDensity::RegularizeDistributionRecursive
(const TUnfoldBinning *binning,ERegMode regmode,
 EDensityMode densityMode,const char *distribution,const char *axisSteering) {
   // recursively regularize distribution(s) using the given settings
   //     binning: distributions for this node an its children are considered
   //     regmode: basic regularisation mode (see class TUnfold)
   //     densityMode: how to apply bin density corrections
   //              (normalisation to bin withd or user factor)
   //     distribution: name of the distribiution where this regularisation
   //             is applied to (if zero, apply to all)
   //     axisSteering: regularisation steering specific to the axes
   //              (see method RegularizeDistribution())
   if((!distribution)|| !TString(distribution).CompareTo(binning->GetName())) {
      RegularizeOneDistribution(binning,regmode,densityMode,axisSteering);
   }
   for(const TUnfoldBinning *child=binning->GetChildNode();child;
       child=child->GetNextNode()) {
      RegularizeDistributionRecursive(child,regmode,densityMode,distribution,
                                      axisSteering);
   }
}

void TUnfoldDensity::RegularizeOneDistribution
(const TUnfoldBinning *binning,ERegMode regmode,
 EDensityMode densityMode,const char *axisSteering)
{
   // regularize the distribution in this node
   //     binning: the distributions to regularize
   //     regmode: basic regularisation mode (see class TUnfold)
   //     densityMode: how to apply bin density corrections
   //              (normalisation to bin withd or user factor)
   //     axisSteering: regularisation steering specific to the axes
   //              (see method RegularizeDistribution()) 
   if(!fRegularisationConditions)
      fRegularisationConditions=new TUnfoldBinning("regularisation");

   TUnfoldBinning *thisRegularisationBinning=
      fRegularisationConditions->AddBinning(binning->GetName());

   // decode steering
   Int_t isOptionGiven[6] = {0};
   binning->DecodeAxisSteering(axisSteering,"uUoObB",isOptionGiven);
   // U implies u
   isOptionGiven[0] |= isOptionGiven[1];
   // O implies o
   isOptionGiven[2] |= isOptionGiven[3];
   // B implies b
   isOptionGiven[4] |= isOptionGiven[5];
#ifdef DEBUG
   cout<<" "<<isOptionGiven[0]
       <<" "<<isOptionGiven[1]
       <<" "<<isOptionGiven[2]
       <<" "<<isOptionGiven[3]
       <<" "<<isOptionGiven[4]
       <<" "<<isOptionGiven[5]
       <<"\n";
#endif
   Info("RegularizeOneDistribution","regularizing %s regMode=%d"
        " densityMode=%d axisSteering=%s",
        binning->GetName(),(Int_t) regmode,(Int_t)densityMode,
        axisSteering ? axisSteering : "");
   Int_t startBin=binning->GetStartBin();
   Int_t endBin=startBin+ binning->GetDistributionNumberOfBins();
   std::vector<Double_t> factor(endBin-startBin);
   Int_t nbin=0;
   for(Int_t bin=startBin;bin<endBin;bin++) {
      factor[bin-startBin]=GetDensityFactor(densityMode,bin);
      if(factor[bin-startBin] !=0.0) nbin++;
   }
#ifdef DEBUG
   cout<<"initial number of bins "<<nbin<<"\n";
#endif
   Int_t dimension=binning->GetDistributionDimension();

   // decide whether to skip underflow/overflow bins
   nbin=0;
   for(Int_t bin=startBin;bin<endBin;bin++) {
      Int_t uStatus,oStatus;
      binning->GetBinUnderflowOverflowStatus(bin,&uStatus,&oStatus);
      if(uStatus & isOptionGiven[1]) factor[bin-startBin]=0.;
      if(oStatus & isOptionGiven[3]) factor[bin-startBin]=0.;
      if(factor[bin-startBin] !=0.0) nbin++;
   }
#ifdef DEBUG
   cout<<"after underflow/overflow bin removal "<<nbin<<"\n";
#endif   
   if(regmode==kRegModeSize) {
      Int_t nRegBins=0;
      // regularize all bins of the distribution, possibly excluding
      // underflow/overflow bins
      for(Int_t bin=startBin;bin<endBin;bin++) {
         if(factor[bin-startBin]==0.0) continue;
         if(AddRegularisationCondition(bin,factor[bin-startBin])) {
            nRegBins++;
         }
      }
      if(nRegBins) {
         thisRegularisationBinning->AddBinning("size",nRegBins);
      }
   } else if((regmode==kRegModeDerivative)||(regmode==kRegModeCurvature)) {
      for(Int_t direction=0;direction<dimension;direction++) {
         // for each direction
         Int_t nRegBins=0;
         Int_t directionMask=(1<<direction);
         Double_t binDistanceNormalisation=
            (isOptionGiven[5] & directionMask)  ?
            binning->GetDistributionAverageBinSize
            (direction,isOptionGiven[0] & directionMask,
             isOptionGiven[2] & directionMask) : 1.0;
         for(Int_t bin=startBin;bin<endBin;bin++) {
            // check whether bin is excluded
            if(factor[bin-startBin]==0.0) continue;
            // for each bin, find the neighbour bins
            Int_t iPrev,iNext;
            Double_t distPrev,distNext;
            binning->GetBinNeighbours
               (bin,direction,&iPrev,&distPrev,&iNext,&distNext);
            if((regmode==kRegModeDerivative)&&(iNext>=0)) {
               Double_t f0 = -factor[bin-startBin];
               Double_t f1 = factor[iNext-startBin];
               if(isOptionGiven[4] & directionMask) {
                  if(distNext>0.0) {
                     f0 *= binDistanceNormalisation/distNext;
                     f1 *= binDistanceNormalisation/distNext;
                  } else {
                     f0=0.;
                     f1=0.;
                  }
               }
               if((f0==0.0)||(f1==0.0)) continue;
               if(AddRegularisationCondition(bin,f0,iNext,f1)) {
                  nRegBins++;
#ifdef DEBUG
                  std::cout<<"Added Reg: bin "<<bin<<" "<<f0
                           <<" next: "<<iNext<<" "<<f1<<"\n";
#endif
               }
            } else if((regmode==kRegModeCurvature)&&(iPrev>=0)&&(iNext>=0)) {
               Double_t f0 = factor[iPrev-startBin];
               Double_t f1 = -factor[bin-startBin];
               Double_t f2 = factor[iNext-startBin];
               if(isOptionGiven[4] & directionMask) {
                  if((distPrev<0.)&&(distNext>0.)) {
                     distPrev= -distPrev;
                     Double_t f=TMath::Power(binDistanceNormalisation,2.)/
                        (distPrev+distNext);
                     f0 *= f/distPrev;
                     f1 *= f*(1./distPrev+1./distNext);
                     f2 *= f/distNext;
                  } else {
                     f0=0.;
                     f1=0.;
                     f2=0.;
                  }
               }
               if((f0==0.0)||(f1==0.0)||(f2==0.0)) continue;
               if(AddRegularisationCondition(iPrev,f0,bin,f1,iNext,f2)) {
                  nRegBins++;
#ifdef DEBUG
                  std::cout<<"Added Reg: prev "<<iPrev<<" "<<f0
                           <<" bin: "<<bin<<" "<<f1
                           <<" next: "<<iNext<<" "<<f2<<"\n";
#endif
               }
            }
         }
         if(nRegBins) {
            TString name;
            if(regmode==kRegModeDerivative) {
               name="derivative_";
            } else if(regmode==kRegModeCurvature) {
               name="curvature_";
            }
            name +=  binning->GetDistributionAxisLabel(direction);
            thisRegularisationBinning->AddBinning(name,nRegBins);
         }
      }
   }
#ifdef DEBUG
   //fLsquared->Print();
#endif
}

TH1 *TUnfoldDensity::GetOutput
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning) const
{
   // retreive unfolding result as histogram
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
     TUnfoldSys::GetOutput(r,binMap);
   }
   if(binMap) {
     delete [] binMap;
   }
   return r;
}

TH1 *TUnfoldDensity::GetBias
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning) const
{
   // retreive unfolding bias as histogram
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetBias(r,binMap);
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetFoldedOutput
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning,Bool_t addBgr) const
{
   // retreive unfolding result folded back by the matrix
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   addBgr: true if the background shall be included
   TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetFoldedOutput(r,binMap);
      if(addBgr) {
         TUnfoldSys::GetBackground(r,0,binMap,0,kFALSE);
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetBackground
(const char *histogramName,const char *bgrSource,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning,
 Int_t includeError,Bool_t clearHist) const
{
   // retreive a background source
   //   histogramName:  name of the histogram
   //   bgrSource: name of the background source
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   include error: +1 if uncorrelated bgr errors should be included
   //                  +2 if correlated bgr errors should be included
   //   clearHist: whether the histogram should be cleared
   //              if false, the background sources are added to the histogram
   TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetBackground(r,bgrSource,binMap,includeError,clearHist);
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetInput
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning) const
{
   // retreive input distribution
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   TUnfoldBinning const *binning=fConstInputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetInput(r,binMap);
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetRhoItotal
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning,TH2 **ematInv) {
   // retreive global correlation coefficients, total error
   // and inverse of error matrix
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   //   ematInv: retreive inverse of error matrix
   //              if ematInv==0 the inverse is not returned
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TH2 *invEmat=0;
      if(ematInv) {
         if(r->GetDimension()==1) {
            TString ematName(histogramName);
            ematName += "_inverseEMAT";
            Int_t *binMap2D=0;
            invEmat=binning->CreateErrorMatrixHistogram
               (ematName,useAxisBinning,&binMap2D,histogramTitle,
                axisSteering);
            if(binMap2D) delete [] binMap2D;
         } else {
            Error("GetRhoItotal",
                  "can not return inverse of error matrix for this binning");
         }
      }
      TUnfoldSys::GetRhoItotal(r,binMap,invEmat);
      if(invEmat) {
         *ematInv=invEmat;
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetRhoIstatbgr
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning,TH2 **ematInv) {
   // retreive global correlation coefficients, input error
   // and inverse of corresponding error matrix
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   //   ematInv: retreive inverse of error matrix
   //              if ematInv==0 the inverse is not returned
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TH2 *invEmat=0;
      if(ematInv) {
         if(r->GetDimension()==1) {
            TString ematName(histogramName);
            ematName += "_inverseEMAT";
            Int_t *binMap2D=0;
            invEmat=binning->CreateErrorMatrixHistogram
               (ematName,useAxisBinning,&binMap2D,histogramTitle,
                axisSteering);
            if(binMap2D) delete [] binMap2D;
         } else {
            Error("GetRhoItotal",
                  "can not return inverse of error matrix for this binning");
         }
      }
      TUnfoldSys::GetRhoI(r,binMap,invEmat);
      if(invEmat) {
         *ematInv=invEmat;
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetDeltaSysSource
(const char *source,const char *histogramName,
 const char *histogramTitle,const char *distributionName,
 const char *axisSteering,Bool_t useAxisBinning) {
   // retreive histogram of systematic 1-sigma shifts
   //   source: name of systematic error
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      if(!TUnfoldSys::GetDeltaSysSource(r,source,binMap)) {
         delete r;
         r=0;
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetDeltaSysBackgroundScale
(const char *bgrSource,const char *histogramName,
 const char *histogramTitle,const char *distributionName,
 const char *axisSteering,Bool_t useAxisBinning) {
   // retreive histogram of systematic 1-sigma shifts due to a background 
   // normalisation uncertainty
   //   source: name of background source
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      if(!TUnfoldSys::GetDeltaSysBackgroundScale(r,bgrSource,binMap)) {
         delete r;
         r=0;
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH1 *TUnfoldDensity::GetDeltaSysTau
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,Bool_t useAxisBinning)
{
   // retreive histogram of systematic 1-sigma shifts due to tau uncertainty
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH1 *r=binning->CreateHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      if(!TUnfoldSys::GetDeltaSysTau(r,binMap)) {
         delete r;
         r=0;
      }
   }
   if(binMap) delete [] binMap;
   return r;
}

TH2 *TUnfoldDensity::GetRhoIJtotal
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning)
{
   // retreive histogram of total corelation coefficients, including systematic
   // uncertainties
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TH2 *r=GetEmatrixTotal
      (histogramName,histogramTitle,distributionName,
       axisSteering,useAxisBinning);
   if(r) {
      for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
         Double_t e_i=r->GetBinContent(i,i);
         if(e_i>0.0) e_i=TMath::Sqrt(e_i);
         else e_i=0.0;
         for(Int_t j=0;j<=r->GetNbinsY()+1;j++) {
            if(i==j) continue;
            Double_t e_j=r->GetBinContent(j,j);
            if(e_j>0.0) e_j=TMath::Sqrt(e_j);
            else e_j=0.0;
            Double_t e_ij=r->GetBinContent(i,j);
            if((e_i>0.0)&&(e_j>0.0)) {
               r->SetBinContent(i,j,e_ij/e_i/e_j);
            } else {
               r->SetBinContent(i,j,0.0);
            }
         }
      }
      for(Int_t i=0;i<=r->GetNbinsX()+1;i++) {
         if(r->GetBinContent(i,i)>0.0) {
            r->SetBinContent(i,i,1.0);
         } else {
            r->SetBinContent(i,i,0.0);
         }
      }
   }
   return r;
}
                           
TH2 *TUnfoldDensity::GetEmatrixSysUncorr
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning)
{
   // get error matrix contribution from uncorrelated errors on the matrix A
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH2 *r=binning->CreateErrorMatrixHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetEmatrixSysUncorr(r,binMap);
   }
   if(binMap) delete [] binMap;
   return r;
}


TH2 *TUnfoldDensity::GetEmatrixSysBackgroundUncorr
(const char *bgrSource,const char *histogramName,
 const char *histogramTitle,const char *distributionName,
 const char *axisSteering,Bool_t useAxisBinning)
{
   // get error matrix from uncorrelated error of one background source
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH2 *r=binning->CreateErrorMatrixHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetEmatrixSysBackgroundUncorr(r,bgrSource,binMap,kFALSE);
   }
   if(binMap) delete [] binMap;
   return r;
}

TH2 *TUnfoldDensity::GetEmatrixInput
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning)
{
   // get error contribution from input vector
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH2 *r=binning->CreateErrorMatrixHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetEmatrixInput(r,binMap);
   }
   if(binMap) delete [] binMap;
   return r;
}

TH2 *TUnfoldDensity::GetProbabilityMatrix
(const char *histogramName,const char *histogramTitle,
 Bool_t useAxisBinning) const
{
   // get matrix of probabilities
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   useAxisBinning: if true, try to get the histogram using
   //                   the original matrix binning
   TH2 *r=TUnfoldBinning::CreateHistogramOfMigrations
      (fConstOutputBins,fConstInputBins,histogramName,
       useAxisBinning,useAxisBinning,histogramTitle);
   TUnfold::GetProbabilityMatrix(r,kHistMapOutputHoriz);
   return r;
}

TH2 *TUnfoldDensity::GetEmatrixTotal
(const char *histogramName,const char *histogramTitle,
 const char *distributionName,const char *axisSteering,
 Bool_t useAxisBinning)
{
   // get total error including systematic,statistical,background,tau errors
   //   histogramName:  name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   distributionName: for complex binning schemes specify the name
   //                of the requested distribution within the TUnfoldBinning
   //                object
   //   axisSteering:
   //       "pattern1;pattern2;...;patternN"
   //       patternI = axis[mode]
   //       axis = name or *
   //       mode = C|U|O 
   //        C: collapse axis into one bin
   //        U: discarde underflow bin
   //        O: discarde overflow bin
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the output histogram
   TUnfoldBinning const *binning=fConstOutputBins->FindNode(distributionName);
   Int_t *binMap=0;
   TH2 *r=binning->CreateErrorMatrixHistogram
      (histogramName,useAxisBinning,&binMap,histogramTitle,axisSteering);
   if(r) {
      TUnfoldSys::GetEmatrixTotal(r,binMap);
   }
   if(binMap) delete [] binMap;
   return r;
}   

TH2 *TUnfoldDensity::GetL
(const char *histogramName,const char *histogramTitle,Bool_t useAxisBinning)
{
   // return the matrix of regularisation conditions in a histogram
   // input:
   //   histogramName: name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   //   useAxisBinning: if true, try to use the axis bin widths
   //                   on the x-axis of the output histogram
   if(fRegularisationConditions && 
      (fRegularisationConditions->GetEndBin()-
       fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
      Warning("GetL",
              "remove invalid scheme of regularisation conditions %d %d",
              fRegularisationConditions->GetEndBin(),fL->GetNrows());
      delete fRegularisationConditions;
      fRegularisationConditions=0;
   }
   if(!fRegularisationConditions) {
      fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
      Warning("GetL","create flat regularisation conditions scheme");
   }
   TH2 *r=TUnfoldBinning::CreateHistogramOfMigrations
      (fConstOutputBins,fRegularisationConditions,histogramName,
       useAxisBinning,useAxisBinning,histogramTitle);
   TUnfold::GetL(r);
   return r;
}

TH1 *TUnfoldDensity::GetLxMinusBias
(const char *histogramName,const char *histogramTitle)
{
   // get regularisation conditions multiplied by result vector minus bias
   //   L(x-biasScale*biasVector)
   // this is a measure of the level of regulartisation required per
   // regularisation condition
   // if there are (negative or positive) spikes,
   // these regularisation conditions dominate
   // over the other regularisation conditions
   // input
   //   histogramName: name of the histogram
   //   histogramTitle: title of the histogram (could be zero)
   TMatrixD dx(*GetX(), TMatrixD::kMinus, fBiasScale * (*fX0));
   TMatrixDSparse *Ldx=MultiplyMSparseM(fL,&dx);
   if(fRegularisationConditions && 
      (fRegularisationConditions->GetEndBin()-
       fRegularisationConditions->GetStartBin()!= fL->GetNrows())) {
      Warning("GetLxMinusBias",
              "remove invalid scheme of regularisation conditions %d %d",
              fRegularisationConditions->GetEndBin(),fL->GetNrows());
      delete fRegularisationConditions;
      fRegularisationConditions=0;
   }
   if(!fRegularisationConditions) {
      fRegularisationConditions=new TUnfoldBinning("regularisation",fL->GetNrows());
      Warning("GetLxMinusBias","create flat regularisation conditions scheme");
   }
   TH1 *r=fRegularisationConditions->CreateHistogram
      (histogramName,kFALSE,0,histogramTitle);
   const Int_t *Ldx_rows=Ldx->GetRowIndexArray();
   const Double_t *Ldx_data=Ldx->GetMatrixArray();
   for(Int_t row=0;row<Ldx->GetNrows();row++) {
      if(Ldx_rows[row]<Ldx_rows[row+1]) {
         r->SetBinContent(row+1,Ldx_data[Ldx_rows[row]]);
      }
   }
   delete Ldx;
   return r;
}

const TUnfoldBinning *TUnfoldDensity::GetInputBinning
(const char *distributionName) const
{
   // find binning scheme, input bins
   //   distributionName : the distribution to locate
   return fConstInputBins->FindNode(distributionName);
}

const TUnfoldBinning *TUnfoldDensity::GetOutputBinning
(const char *distributionName) const
{
   // find binning scheme, output bins
   //   distributionName : the distribution to locate
   return fConstOutputBins->FindNode(distributionName);
}

Int_t TUnfoldDensity::ScanTau
(Int_t nPoint,Double_t tauMin,Double_t tauMax,TSpline **scanResult,
 Int_t mode,const char *distribution,const char *axisSteering,
 TGraph **lCurvePlot,TSpline **logTauXPlot,TSpline **logTauYPlot)
{
   // scan some variable as a function of tau and determine the minimum
   // input:
   //   nPoint: number of points to be scanned on the resulting curve
   //   tauMin: smallest tau value to study
   //   tauMax: largest tau value to study
   //     if (mauMin,tauMax) do not correspond to a valid tau range
   //     (e.g. tauMin=tauMax=0.0) then the tau range is determined
   //     automatically
   //   mode,distribution,axisSteering: argument to GetScanVariable()
   // output:
   //   scanResult: output spline of the variable as a function of tau
   // the following plots are produced on request (if pointers are non-zero)
   //   lCurvePlot: for monitoring: the L-curve
   //   logTauXPlot: for monitoring: L-curve(x) as a function of log(tau)
   //   logTauYPlot: for monitoring: L-curve(y) as a function of log(tau)
   // return value: the coordinate number (0..nPoint-1) corresponding to the
   //   final choice of tau
   typedef std::map<Double_t,Double_t> TauScan_t;
   typedef std::map<Double_t,std::pair<Double_t,Double_t> > LCurve_t;
   TauScan_t curve;
   LCurve_t lcurve;

   //==========================================================
   // algorithm:
   //  (1) do the unfolding for nPoint-1 points
   //      and store the results in the map
   //        curve
   //    (1a) store minimum and maximum tau to curve
   //    (1b) insert additional points, until nPoint-1 values
   //          have been calculated
   //
   //  (2) determine the best choice of tau
   //      do the unfolding for this point
   //      and store the result in
   //        curve
   //  (3) return the result in scanResult
   
   //==========================================================
   //  (1) do the unfolding for nPoint-1 points
   //      and store the results in
   //        curve
   //    (1a) store minimum and maximum tau to curve
   
   if((tauMin<=0)||(tauMax<=0.0)||(tauMin>=tauMax)) {
      // here no range is given, has to be determined automatically
      // the maximum tau is determined from the chi**2 values 
      // observed from unfolding without regulatisation
      
      // first unfolding, without regularisation
      DoUnfold(0.0);

      // if the number of degrees of freedom is too small, create an error
      if(GetNdf()<=0) {
         Error("ScanTau","too few input bins, NDF<=0 %d",GetNdf());  
      }
      Double_t X0=GetLcurveX();
      Double_t Y0=GetLcurveY();
      Double_t y0=GetScanVariable(mode,distribution,axisSteering);
      Info("ScanTau","logtau=-Infinity y=%lf X=%lf Y=%lf",y0,X0,Y0);
      {
         // unfolding guess maximum tau and store it
         Double_t logTau=
            0.5*(TMath::Log10(GetChi2A()+3.*TMath::Sqrt(GetNdf()+1.0))
                 -GetLcurveY());
         DoUnfold(TMath::Power(10.,logTau));
         if((!TMath::Finite(GetLcurveX())) ||(!TMath::Finite(GetLcurveY()))) {
            Fatal("ScanTau","problem (missing regularisation?) X=%f Y=%f",
                  GetLcurveX(),GetLcurveY());
         }
         Double_t y=GetScanVariable(mode,distribution,axisSteering);
         curve[logTau]=y;
         lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
         Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
              GetLcurveX(),GetLcurveY());
      }
      // minimum tau is chosen such that it is less than
      // 1% different from the case of no regularisation
      // here, several points are inserted as needed
      while(((int)curve.size()<nPoint-1)&&
            ((TMath::Abs(GetLcurveX()-X0)>0.00432)||
             (TMath::Abs(GetLcurveY()-Y0)>0.00432))) {
         Double_t logTau=(*curve.begin()).first-0.5;
         DoUnfold(TMath::Power(10.,logTau));
         Double_t y=GetScanVariable(mode,distribution,axisSteering);
         curve[logTau]=y;
         lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
         Info("ScanTay","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
              GetLcurveX(),GetLcurveY());
      }
   } else {
      Double_t logTauMin=TMath::Log10(tauMin);
      Double_t logTauMax=TMath::Log10(tauMax);
      if(nPoint>1) {
         // insert maximum tau
         DoUnfold(TMath::Power(10.,logTauMax));
         Double_t y=GetScanVariable(mode,distribution,axisSteering);
         curve[logTauMax]=y;
         lcurve[logTauMax]=std::make_pair(GetLcurveX(),GetLcurveY());
         Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMax,y,
              GetLcurveX(),GetLcurveY());
      }
      // insert minimum tau
      DoUnfold(TMath::Power(10.,logTauMin));
      Double_t y=GetScanVariable(mode,distribution,axisSteering);
      curve[logTauMin]=y;
      lcurve[logTauMin]=std::make_pair(GetLcurveX(),GetLcurveY());
      Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTauMin,y,
           GetLcurveX(),GetLcurveY());
   }

   //==========================================================
   //    (1b) insert additional points, until nPoint-1 values
   //          have been calculated
   while((int)curve.size()<nPoint-1) {
      // insert additional points
      // points are inserted such that the largest interval in log(tau)
      // is divided into two smaller intervals
      // however, there is a penalty term for subdividing intervals
      // which are far away from the minimum
      TauScan_t::const_iterator i0,i1;
      i0=curve.begin();
      // locate minimum
      Double_t logTauYMin=(*i0).first;
      Double_t yMin=(*i0).second;
      for(;i0!=curve.end();i0++) {
         if((*i0).second<yMin) {
            yMin=(*i0).second;
            logTauYMin=(*i0).first;
         }
      }
      // insert additional points such that large log(tau) intervals
      // near the minimum rho are divided into two
      i0=curve.begin();
      i1=i0;
      Double_t distMax=0.0;
      Double_t logTau=0.0;
      for(i1++;i1!=curve.end();i1++) {
         Double_t dist;
         // check size of rho interval
         dist=TMath::Abs((*i0).first-(*i1).first)
            // penalty term if distance from rhoMax is large
            +0.25*TMath::Power(0.5*((*i0).first+(*i1).first)-logTauYMin,2.)/
            ((*curve.rbegin()).first-(*curve.begin()).first)/nPoint;
         if((dist<=0.0)||(dist>distMax)) {
            distMax=dist;
            logTau=0.5*((*i0).first+(*i1).first);
         }
         i0=i1;
      }
      DoUnfold(TMath::Power(10.,logTau));
      Double_t y=GetScanVariable(mode,distribution,axisSteering);
      curve[logTau]=y;
      lcurve[logTau]=std::make_pair(GetLcurveX(),GetLcurveY());
      Info("ScanTau","logtau=%lf y=%lf X=%lf Y=%lf",logTau,y,
           GetLcurveX(),GetLcurveY());
   }

   //==========================================================
   //  (2) determine the best choice of tau
   //      do the unfolding for this point
   //      and store the result in
   //        curve

   Double_t cTmin=0.0;
   {
   Double_t *cTi=new Double_t[curve.size()];
   Double_t *cCi=new Double_t[curve.size()];
   Int_t n=0;
   for(TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
      cTi[n]=(*i).first;
      cCi[n]=(*i).second;
      n++;
   }
   // create rho Spline
   TSpline3 *splineC=new TSpline3("L curve curvature",cTi,cCi,n);
   // find the maximum of the curvature
   // if the parameter iskip is non-zero, then iskip points are
   // ignored when looking for the largest curvature
   // (there are problems with the curvature determined from the first
   //  few points of splineX,splineY in the algorithm above)
   Int_t iskip=0;
   if(n>3) iskip=1;
   if(n>6) iskip=2;
   Double_t cCmin=cCi[iskip];
   cTmin=cTi[iskip];
   for(Int_t i=iskip;i<n-1-iskip;i++) {
      // find minimum on this spline section
      // check boundary conditions for x[i+1]
      Double_t xMin=cTi[i+1];
      Double_t yMin=cCi[i+1];
      if(cCi[i]<yMin) {
         yMin=cCi[i];
         xMin=cTi[i];
      }
      // find minimum for x[i]<x<x[i+1]
      // get spline coefficiencts and solve equation
      //   derivative(x)==0
      Double_t x,y,b,c,d;
      splineC->GetCoeff(i,x,y,b,c,d);
      // coefficiencts of quadratic equation
      Double_t m_p_half=-c/(3.*d);
      Double_t q=b/(3.*d);
      Double_t discr=m_p_half*m_p_half-q;
      if(discr>=0.0) {
         // solution found
         discr=TMath::Sqrt(discr);
         Double_t xx;
         if(m_p_half>0.0) {
            xx = m_p_half + discr;
         } else {
            xx = m_p_half - discr;
         }
         Double_t dx=cTi[i+1]-x;
         // check first solution
         if((xx>0.0)&&(xx<dx)) {
            y=splineC->Eval(x+xx);
            if(y<yMin) {
               yMin=y;
               xMin=x+xx;
            }
         }
         // second solution
         if(xx !=0.0) {
            xx= q/xx;
         } else {
            xx=0.0;
         }
         // check second solution
         if((xx>0.0)&&(xx<dx)) {
            y=splineC->Eval(x+xx);
            if(y<yMin) {
               yMin=y;
               xMin=x+xx;
            }
         }
      }
      // check whether this local minimum is a global minimum
      if(yMin<cCmin) {
         cCmin=yMin;
         cTmin=xMin;
      }
   }
   delete splineC;
   delete[] cTi;
   delete[] cCi;
   }
   Double_t logTauFin=cTmin;
   DoUnfold(TMath::Power(10.,logTauFin));
   {
      Double_t y=GetScanVariable(mode,distribution,axisSteering);
      curve[logTauFin]=y;
      lcurve[logTauFin]=std::make_pair(GetLcurveX(),GetLcurveY());
      Info("ScanTau","Result logtau=%lf y=%lf X=%lf Y=%lf",logTauFin,y,
           GetLcurveX(),GetLcurveY());
   }
   //==========================================================
   //  (3) return the result in
   //       scanResult lCurve logTauX logTauY

   Int_t bestChoice=-1;
   if(curve.size()>0) {
      Double_t *y=new Double_t[curve.size()];
      Double_t *logT=new Double_t[curve.size()];
      int n=0;
      for( TauScan_t::const_iterator i=curve.begin();i!=curve.end();i++) {
         if(logTauFin==(*i).first) {
            bestChoice=n;
         }
         y[n]=(*i).second;
         logT[n]=(*i).first;
         n++;
      }
      if(scanResult) {
         TString name;
         name = TString::Format("scan(%d,",mode);
         if(distribution) name+= distribution;
         name += ",";
         if(axisSteering) name += axisSteering;
         name +=")";
         (*scanResult)=new TSpline3(name+"%log(tau)",logT,y,n);
      }
      delete[] y;
      delete[] logT;
   }
   if(lcurve.size()) {
      Double_t *logT=new Double_t[lcurve.size()];
      Double_t *x=new Double_t[lcurve.size()];
      Double_t *y=new Double_t[lcurve.size()];
      Int_t n=0;
      for(LCurve_t::const_iterator i=lcurve.begin();i!=lcurve.end();i++) {
         logT[n]=(*i).first;
         x[n]=(*i).second.first;
         y[n]=(*i).second.second;
         //cout<<logT[n]<<" "<< x[n]<<" "<<y[n]<<"\n";
         n++;
      }
      if(lCurvePlot) {
         *lCurvePlot=new TGraph(n,x,y); 
         (*lCurvePlot)->SetTitle("L curve");
      }
      if(logTauXPlot)
         *logTauXPlot=new TSpline3("log(chi**2)%log(tau)",logT,x,n);
      if(logTauYPlot)
         *logTauYPlot=new TSpline3("log(reg.cond)%log(tau)",logT,y,n);
      delete [] y;
      delete [] x;
      delete [] logT;
   }
   return bestChoice;
}

Double_t TUnfoldDensity::GetScanVariable
(Int_t mode,const char *distribution,const char *axisSteering)
{
   // calculate variable for ScanTau()
   // the unfolding is repeated for various choices of tau.
   // For each tau, after unfolding, the ScanTau() method calls
   // GetScanVariable() to determine the value of the variable which
   // is to be scanned
   //
   // the variable is expected to have a minimum near the "optimal" choice
   // of tau
   //
   // input:
   //    mode : define the type of variable to be calculated
   //    distribution : define the distribution for which the variable
   //              is to be calculated
   //        the following modes are implemented:
   //          kEScanTauRhoAvg : average global correlation from input data
   //          kEScanTauRhoSquaredAvg : average global correlation squared
   //                                   from input data
   //          kEScanTauRhoMax : maximum global correlation from input data
   //          kEScanTauRhoAvgSys : average global correlation
   //                                 including systematic uncertainties
   //          kEScanTauRhoAvgSquaredSys : average global correlation squared
   //                                 including systematic uncertainties
   //          kEScanTauRhoMaxSys : maximum global correlation
   //                                 including systematic uncertainties
   //    distribution : name of the TUnfoldBinning node
   //                   for which to calculate the correlations
   //    axisSteering : axis steering for calculating the correlations
   //              the distribution
   // return: the value of the variable as determined from the present
   //    unfolding

   Double_t r=0.0;
   TString name("GetScanVariable(");
   name += TString::Format("%d,",mode);
   if(distribution) name += distribution;
   name += ",";
   if(axisSteering) name += axisSteering;
   name += ")";
   TH1 *rhoi=0;
   if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoMax)||
      (mode==kEScanTauRhoSquareAvg)) {
      rhoi=GetRhoIstatbgr(name,0,distribution,axisSteering,kFALSE);
   } else if((mode==kEScanTauRhoAvgSys)||(mode==kEScanTauRhoMaxSys)||
             (mode==kEScanTauRhoSquareAvgSys)) {
      rhoi=GetRhoItotal(name,0,distribution,axisSteering,kFALSE);
   }
   if(rhoi) {
      Double_t sum=0.0;
      Double_t sumSquare=0.0;
      Double_t rhoMax=0.0;
      Int_t n=0;
      for(Int_t i=0;i<=rhoi->GetNbinsX()+1;i++) {
         Double_t c=rhoi->GetBinContent(i);
         if(c>=0.) {
            if(c>rhoMax) rhoMax=c;
            sum += c;
            sumSquare += c*c;
            n ++;
         }
      }
      if((mode==kEScanTauRhoAvg)||(mode==kEScanTauRhoAvgSys)) {
         r=sum/n;
      } else if((mode==kEScanTauRhoSquareAvg)||
                (mode==kEScanTauRhoSquareAvgSys)) {
         r=sum/n;
      } else {
         r=rhoMax;
      }
      // cout<<r<<" "<<GetRhoAvg()<<" "<<GetRhoMax()<<"\n";
      delete rhoi;
   } else {
      Fatal("GetScanVariable","mode %d not implemented",mode);
   }
   return r;
}