// 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;
}
 TUnfoldDensity.cxx:1
 TUnfoldDensity.cxx:2
 TUnfoldDensity.cxx:3
 TUnfoldDensity.cxx:4
 TUnfoldDensity.cxx:5
 TUnfoldDensity.cxx:6
 TUnfoldDensity.cxx:7
 TUnfoldDensity.cxx:8
 TUnfoldDensity.cxx:9
 TUnfoldDensity.cxx:10
 TUnfoldDensity.cxx:11
 TUnfoldDensity.cxx:12
 TUnfoldDensity.cxx:13
 TUnfoldDensity.cxx:14
 TUnfoldDensity.cxx:15
 TUnfoldDensity.cxx:16
 TUnfoldDensity.cxx:17
 TUnfoldDensity.cxx:18
 TUnfoldDensity.cxx:19
 TUnfoldDensity.cxx:20
 TUnfoldDensity.cxx:21
 TUnfoldDensity.cxx:22
 TUnfoldDensity.cxx:23
 TUnfoldDensity.cxx:24
 TUnfoldDensity.cxx:25
 TUnfoldDensity.cxx:26
 TUnfoldDensity.cxx:27
 TUnfoldDensity.cxx:28
 TUnfoldDensity.cxx:29
 TUnfoldDensity.cxx:30
 TUnfoldDensity.cxx:31
 TUnfoldDensity.cxx:32
 TUnfoldDensity.cxx:33
 TUnfoldDensity.cxx:34
 TUnfoldDensity.cxx:35
 TUnfoldDensity.cxx:36
 TUnfoldDensity.cxx:37
 TUnfoldDensity.cxx:38
 TUnfoldDensity.cxx:39
 TUnfoldDensity.cxx:40
 TUnfoldDensity.cxx:41
 TUnfoldDensity.cxx:42
 TUnfoldDensity.cxx:43
 TUnfoldDensity.cxx:44
 TUnfoldDensity.cxx:45
 TUnfoldDensity.cxx:46
 TUnfoldDensity.cxx:47
 TUnfoldDensity.cxx:48
 TUnfoldDensity.cxx:49
 TUnfoldDensity.cxx:50
 TUnfoldDensity.cxx:51
 TUnfoldDensity.cxx:52
 TUnfoldDensity.cxx:53
 TUnfoldDensity.cxx:54
 TUnfoldDensity.cxx:55
 TUnfoldDensity.cxx:56
 TUnfoldDensity.cxx:57
 TUnfoldDensity.cxx:58
 TUnfoldDensity.cxx:59
 TUnfoldDensity.cxx:60
 TUnfoldDensity.cxx:61
 TUnfoldDensity.cxx:62
 TUnfoldDensity.cxx:63
 TUnfoldDensity.cxx:64
 TUnfoldDensity.cxx:65
 TUnfoldDensity.cxx:66
 TUnfoldDensity.cxx:67
 TUnfoldDensity.cxx:68
 TUnfoldDensity.cxx:69
 TUnfoldDensity.cxx:70
 TUnfoldDensity.cxx:71
 TUnfoldDensity.cxx:72
 TUnfoldDensity.cxx:73
 TUnfoldDensity.cxx:74
 TUnfoldDensity.cxx:75
 TUnfoldDensity.cxx:76
 TUnfoldDensity.cxx:77
 TUnfoldDensity.cxx:78
 TUnfoldDensity.cxx:79
 TUnfoldDensity.cxx:80
 TUnfoldDensity.cxx:81
 TUnfoldDensity.cxx:82
 TUnfoldDensity.cxx:83
 TUnfoldDensity.cxx:84
 TUnfoldDensity.cxx:85
 TUnfoldDensity.cxx:86
 TUnfoldDensity.cxx:87
 TUnfoldDensity.cxx:88
 TUnfoldDensity.cxx:89
 TUnfoldDensity.cxx:90
 TUnfoldDensity.cxx:91
 TUnfoldDensity.cxx:92
 TUnfoldDensity.cxx:93
 TUnfoldDensity.cxx:94
 TUnfoldDensity.cxx:95
 TUnfoldDensity.cxx:96
 TUnfoldDensity.cxx:97
 TUnfoldDensity.cxx:98
 TUnfoldDensity.cxx:99
 TUnfoldDensity.cxx:100
 TUnfoldDensity.cxx:101
 TUnfoldDensity.cxx:102
 TUnfoldDensity.cxx:103
 TUnfoldDensity.cxx:104
 TUnfoldDensity.cxx:105
 TUnfoldDensity.cxx:106
 TUnfoldDensity.cxx:107
 TUnfoldDensity.cxx:108
 TUnfoldDensity.cxx:109
 TUnfoldDensity.cxx:110
 TUnfoldDensity.cxx:111
 TUnfoldDensity.cxx:112
 TUnfoldDensity.cxx:113
 TUnfoldDensity.cxx:114
 TUnfoldDensity.cxx:115
 TUnfoldDensity.cxx:116
 TUnfoldDensity.cxx:117
 TUnfoldDensity.cxx:118
 TUnfoldDensity.cxx:119
 TUnfoldDensity.cxx:120
 TUnfoldDensity.cxx:121
 TUnfoldDensity.cxx:122
 TUnfoldDensity.cxx:123
 TUnfoldDensity.cxx:124
 TUnfoldDensity.cxx:125
 TUnfoldDensity.cxx:126
 TUnfoldDensity.cxx:127
 TUnfoldDensity.cxx:128
 TUnfoldDensity.cxx:129
 TUnfoldDensity.cxx:130
 TUnfoldDensity.cxx:131
 TUnfoldDensity.cxx:132
 TUnfoldDensity.cxx:133
 TUnfoldDensity.cxx:134
 TUnfoldDensity.cxx:135
 TUnfoldDensity.cxx:136
 TUnfoldDensity.cxx:137
 TUnfoldDensity.cxx:138
 TUnfoldDensity.cxx:139
 TUnfoldDensity.cxx:140
 TUnfoldDensity.cxx:141
 TUnfoldDensity.cxx:142
 TUnfoldDensity.cxx:143
 TUnfoldDensity.cxx:144
 TUnfoldDensity.cxx:145
 TUnfoldDensity.cxx:146
 TUnfoldDensity.cxx:147
 TUnfoldDensity.cxx:148
 TUnfoldDensity.cxx:149
 TUnfoldDensity.cxx:150
 TUnfoldDensity.cxx:151
 TUnfoldDensity.cxx:152
 TUnfoldDensity.cxx:153
 TUnfoldDensity.cxx:154
 TUnfoldDensity.cxx:155
 TUnfoldDensity.cxx:156
 TUnfoldDensity.cxx:157
 TUnfoldDensity.cxx:158
 TUnfoldDensity.cxx:159
 TUnfoldDensity.cxx:160
 TUnfoldDensity.cxx:161
 TUnfoldDensity.cxx:162
 TUnfoldDensity.cxx:163
 TUnfoldDensity.cxx:164
 TUnfoldDensity.cxx:165
 TUnfoldDensity.cxx:166
 TUnfoldDensity.cxx:167
 TUnfoldDensity.cxx:168
 TUnfoldDensity.cxx:169
 TUnfoldDensity.cxx:170
 TUnfoldDensity.cxx:171
 TUnfoldDensity.cxx:172
 TUnfoldDensity.cxx:173
 TUnfoldDensity.cxx:174
 TUnfoldDensity.cxx:175
 TUnfoldDensity.cxx:176
 TUnfoldDensity.cxx:177
 TUnfoldDensity.cxx:178
 TUnfoldDensity.cxx:179
 TUnfoldDensity.cxx:180
 TUnfoldDensity.cxx:181
 TUnfoldDensity.cxx:182
 TUnfoldDensity.cxx:183
 TUnfoldDensity.cxx:184
 TUnfoldDensity.cxx:185
 TUnfoldDensity.cxx:186
 TUnfoldDensity.cxx:187
 TUnfoldDensity.cxx:188
 TUnfoldDensity.cxx:189
 TUnfoldDensity.cxx:190
 TUnfoldDensity.cxx:191
 TUnfoldDensity.cxx:192
 TUnfoldDensity.cxx:193
 TUnfoldDensity.cxx:194
 TUnfoldDensity.cxx:195
 TUnfoldDensity.cxx:196
 TUnfoldDensity.cxx:197
 TUnfoldDensity.cxx:198
 TUnfoldDensity.cxx:199
 TUnfoldDensity.cxx:200
 TUnfoldDensity.cxx:201
 TUnfoldDensity.cxx:202
 TUnfoldDensity.cxx:203
 TUnfoldDensity.cxx:204
 TUnfoldDensity.cxx:205
 TUnfoldDensity.cxx:206
 TUnfoldDensity.cxx:207
 TUnfoldDensity.cxx:208
 TUnfoldDensity.cxx:209
 TUnfoldDensity.cxx:210
 TUnfoldDensity.cxx:211
 TUnfoldDensity.cxx:212
 TUnfoldDensity.cxx:213
 TUnfoldDensity.cxx:214
 TUnfoldDensity.cxx:215
 TUnfoldDensity.cxx:216
 TUnfoldDensity.cxx:217
 TUnfoldDensity.cxx:218
 TUnfoldDensity.cxx:219
 TUnfoldDensity.cxx:220
 TUnfoldDensity.cxx:221
 TUnfoldDensity.cxx:222
 TUnfoldDensity.cxx:223
 TUnfoldDensity.cxx:224
 TUnfoldDensity.cxx:225
 TUnfoldDensity.cxx:226
 TUnfoldDensity.cxx:227
 TUnfoldDensity.cxx:228
 TUnfoldDensity.cxx:229
 TUnfoldDensity.cxx:230
 TUnfoldDensity.cxx:231
 TUnfoldDensity.cxx:232
 TUnfoldDensity.cxx:233
 TUnfoldDensity.cxx:234
 TUnfoldDensity.cxx:235
 TUnfoldDensity.cxx:236
 TUnfoldDensity.cxx:237
 TUnfoldDensity.cxx:238
 TUnfoldDensity.cxx:239
 TUnfoldDensity.cxx:240
 TUnfoldDensity.cxx:241
 TUnfoldDensity.cxx:242
 TUnfoldDensity.cxx:243
 TUnfoldDensity.cxx:244
 TUnfoldDensity.cxx:245
 TUnfoldDensity.cxx:246
 TUnfoldDensity.cxx:247
 TUnfoldDensity.cxx:248
 TUnfoldDensity.cxx:249
 TUnfoldDensity.cxx:250
 TUnfoldDensity.cxx:251
 TUnfoldDensity.cxx:252
 TUnfoldDensity.cxx:253
 TUnfoldDensity.cxx:254
 TUnfoldDensity.cxx:255
 TUnfoldDensity.cxx:256
 TUnfoldDensity.cxx:257
 TUnfoldDensity.cxx:258
 TUnfoldDensity.cxx:259
 TUnfoldDensity.cxx:260
 TUnfoldDensity.cxx:261
 TUnfoldDensity.cxx:262
 TUnfoldDensity.cxx:263
 TUnfoldDensity.cxx:264
 TUnfoldDensity.cxx:265
 TUnfoldDensity.cxx:266
 TUnfoldDensity.cxx:267
 TUnfoldDensity.cxx:268
 TUnfoldDensity.cxx:269
 TUnfoldDensity.cxx:270
 TUnfoldDensity.cxx:271
 TUnfoldDensity.cxx:272
 TUnfoldDensity.cxx:273
 TUnfoldDensity.cxx:274
 TUnfoldDensity.cxx:275
 TUnfoldDensity.cxx:276
 TUnfoldDensity.cxx:277
 TUnfoldDensity.cxx:278
 TUnfoldDensity.cxx:279
 TUnfoldDensity.cxx:280
 TUnfoldDensity.cxx:281
 TUnfoldDensity.cxx:282
 TUnfoldDensity.cxx:283
 TUnfoldDensity.cxx:284
 TUnfoldDensity.cxx:285
 TUnfoldDensity.cxx:286
 TUnfoldDensity.cxx:287
 TUnfoldDensity.cxx:288
 TUnfoldDensity.cxx:289
 TUnfoldDensity.cxx:290
 TUnfoldDensity.cxx:291
 TUnfoldDensity.cxx:292
 TUnfoldDensity.cxx:293
 TUnfoldDensity.cxx:294
 TUnfoldDensity.cxx:295
 TUnfoldDensity.cxx:296
 TUnfoldDensity.cxx:297
 TUnfoldDensity.cxx:298
 TUnfoldDensity.cxx:299
 TUnfoldDensity.cxx:300
 TUnfoldDensity.cxx:301
 TUnfoldDensity.cxx:302
 TUnfoldDensity.cxx:303
 TUnfoldDensity.cxx:304
 TUnfoldDensity.cxx:305
 TUnfoldDensity.cxx:306
 TUnfoldDensity.cxx:307
 TUnfoldDensity.cxx:308
 TUnfoldDensity.cxx:309
 TUnfoldDensity.cxx:310
 TUnfoldDensity.cxx:311
 TUnfoldDensity.cxx:312
 TUnfoldDensity.cxx:313
 TUnfoldDensity.cxx:314
 TUnfoldDensity.cxx:315
 TUnfoldDensity.cxx:316
 TUnfoldDensity.cxx:317
 TUnfoldDensity.cxx:318
 TUnfoldDensity.cxx:319
 TUnfoldDensity.cxx:320
 TUnfoldDensity.cxx:321
 TUnfoldDensity.cxx:322
 TUnfoldDensity.cxx:323
 TUnfoldDensity.cxx:324
 TUnfoldDensity.cxx:325
 TUnfoldDensity.cxx:326
 TUnfoldDensity.cxx:327
 TUnfoldDensity.cxx:328
 TUnfoldDensity.cxx:329
 TUnfoldDensity.cxx:330
 TUnfoldDensity.cxx:331
 TUnfoldDensity.cxx:332
 TUnfoldDensity.cxx:333
 TUnfoldDensity.cxx:334
 TUnfoldDensity.cxx:335
 TUnfoldDensity.cxx:336
 TUnfoldDensity.cxx:337
 TUnfoldDensity.cxx:338
 TUnfoldDensity.cxx:339
 TUnfoldDensity.cxx:340
 TUnfoldDensity.cxx:341
 TUnfoldDensity.cxx:342
 TUnfoldDensity.cxx:343
 TUnfoldDensity.cxx:344
 TUnfoldDensity.cxx:345
 TUnfoldDensity.cxx:346
 TUnfoldDensity.cxx:347
 TUnfoldDensity.cxx:348
 TUnfoldDensity.cxx:349
 TUnfoldDensity.cxx:350
 TUnfoldDensity.cxx:351
 TUnfoldDensity.cxx:352
 TUnfoldDensity.cxx:353
 TUnfoldDensity.cxx:354
 TUnfoldDensity.cxx:355
 TUnfoldDensity.cxx:356
 TUnfoldDensity.cxx:357
 TUnfoldDensity.cxx:358
 TUnfoldDensity.cxx:359
 TUnfoldDensity.cxx:360
 TUnfoldDensity.cxx:361
 TUnfoldDensity.cxx:362
 TUnfoldDensity.cxx:363
 TUnfoldDensity.cxx:364
 TUnfoldDensity.cxx:365
 TUnfoldDensity.cxx:366
 TUnfoldDensity.cxx:367
 TUnfoldDensity.cxx:368
 TUnfoldDensity.cxx:369
 TUnfoldDensity.cxx:370
 TUnfoldDensity.cxx:371
 TUnfoldDensity.cxx:372
 TUnfoldDensity.cxx:373
 TUnfoldDensity.cxx:374
 TUnfoldDensity.cxx:375
 TUnfoldDensity.cxx:376
 TUnfoldDensity.cxx:377
 TUnfoldDensity.cxx:378
 TUnfoldDensity.cxx:379
 TUnfoldDensity.cxx:380
 TUnfoldDensity.cxx:381
 TUnfoldDensity.cxx:382
 TUnfoldDensity.cxx:383
 TUnfoldDensity.cxx:384
 TUnfoldDensity.cxx:385
 TUnfoldDensity.cxx:386
 TUnfoldDensity.cxx:387
 TUnfoldDensity.cxx:388
 TUnfoldDensity.cxx:389
 TUnfoldDensity.cxx:390
 TUnfoldDensity.cxx:391
 TUnfoldDensity.cxx:392
 TUnfoldDensity.cxx:393
 TUnfoldDensity.cxx:394
 TUnfoldDensity.cxx:395
 TUnfoldDensity.cxx:396
 TUnfoldDensity.cxx:397
 TUnfoldDensity.cxx:398
 TUnfoldDensity.cxx:399
 TUnfoldDensity.cxx:400
 TUnfoldDensity.cxx:401
 TUnfoldDensity.cxx:402
 TUnfoldDensity.cxx:403
 TUnfoldDensity.cxx:404
 TUnfoldDensity.cxx:405
 TUnfoldDensity.cxx:406
 TUnfoldDensity.cxx:407
 TUnfoldDensity.cxx:408
 TUnfoldDensity.cxx:409
 TUnfoldDensity.cxx:410
 TUnfoldDensity.cxx:411
 TUnfoldDensity.cxx:412
 TUnfoldDensity.cxx:413
 TUnfoldDensity.cxx:414
 TUnfoldDensity.cxx:415
 TUnfoldDensity.cxx:416
 TUnfoldDensity.cxx:417
 TUnfoldDensity.cxx:418
 TUnfoldDensity.cxx:419
 TUnfoldDensity.cxx:420
 TUnfoldDensity.cxx:421
 TUnfoldDensity.cxx:422
 TUnfoldDensity.cxx:423
 TUnfoldDensity.cxx:424
 TUnfoldDensity.cxx:425
 TUnfoldDensity.cxx:426
 TUnfoldDensity.cxx:427
 TUnfoldDensity.cxx:428
 TUnfoldDensity.cxx:429
 TUnfoldDensity.cxx:430
 TUnfoldDensity.cxx:431
 TUnfoldDensity.cxx:432
 TUnfoldDensity.cxx:433
 TUnfoldDensity.cxx:434
 TUnfoldDensity.cxx:435
 TUnfoldDensity.cxx:436
 TUnfoldDensity.cxx:437
 TUnfoldDensity.cxx:438
 TUnfoldDensity.cxx:439
 TUnfoldDensity.cxx:440
 TUnfoldDensity.cxx:441
 TUnfoldDensity.cxx:442
 TUnfoldDensity.cxx:443
 TUnfoldDensity.cxx:444
 TUnfoldDensity.cxx:445
 TUnfoldDensity.cxx:446
 TUnfoldDensity.cxx:447
 TUnfoldDensity.cxx:448
 TUnfoldDensity.cxx:449
 TUnfoldDensity.cxx:450
 TUnfoldDensity.cxx:451
 TUnfoldDensity.cxx:452
 TUnfoldDensity.cxx:453
 TUnfoldDensity.cxx:454
 TUnfoldDensity.cxx:455
 TUnfoldDensity.cxx:456
 TUnfoldDensity.cxx:457
 TUnfoldDensity.cxx:458
 TUnfoldDensity.cxx:459
 TUnfoldDensity.cxx:460
 TUnfoldDensity.cxx:461
 TUnfoldDensity.cxx:462
 TUnfoldDensity.cxx:463
 TUnfoldDensity.cxx:464
 TUnfoldDensity.cxx:465
 TUnfoldDensity.cxx:466
 TUnfoldDensity.cxx:467
 TUnfoldDensity.cxx:468
 TUnfoldDensity.cxx:469
 TUnfoldDensity.cxx:470
 TUnfoldDensity.cxx:471
 TUnfoldDensity.cxx:472
 TUnfoldDensity.cxx:473
 TUnfoldDensity.cxx:474
 TUnfoldDensity.cxx:475
 TUnfoldDensity.cxx:476
 TUnfoldDensity.cxx:477
 TUnfoldDensity.cxx:478
 TUnfoldDensity.cxx:479
 TUnfoldDensity.cxx:480
 TUnfoldDensity.cxx:481
 TUnfoldDensity.cxx:482
 TUnfoldDensity.cxx:483
 TUnfoldDensity.cxx:484
 TUnfoldDensity.cxx:485
 TUnfoldDensity.cxx:486
 TUnfoldDensity.cxx:487
 TUnfoldDensity.cxx:488
 TUnfoldDensity.cxx:489
 TUnfoldDensity.cxx:490
 TUnfoldDensity.cxx:491
 TUnfoldDensity.cxx:492
 TUnfoldDensity.cxx:493
 TUnfoldDensity.cxx:494
 TUnfoldDensity.cxx:495
 TUnfoldDensity.cxx:496
 TUnfoldDensity.cxx:497
 TUnfoldDensity.cxx:498
 TUnfoldDensity.cxx:499
 TUnfoldDensity.cxx:500
 TUnfoldDensity.cxx:501
 TUnfoldDensity.cxx:502
 TUnfoldDensity.cxx:503
 TUnfoldDensity.cxx:504
 TUnfoldDensity.cxx:505
 TUnfoldDensity.cxx:506
 TUnfoldDensity.cxx:507
 TUnfoldDensity.cxx:508
 TUnfoldDensity.cxx:509
 TUnfoldDensity.cxx:510
 TUnfoldDensity.cxx:511
 TUnfoldDensity.cxx:512
 TUnfoldDensity.cxx:513
 TUnfoldDensity.cxx:514
 TUnfoldDensity.cxx:515
 TUnfoldDensity.cxx:516
 TUnfoldDensity.cxx:517
 TUnfoldDensity.cxx:518
 TUnfoldDensity.cxx:519
 TUnfoldDensity.cxx:520
 TUnfoldDensity.cxx:521
 TUnfoldDensity.cxx:522
 TUnfoldDensity.cxx:523
 TUnfoldDensity.cxx:524
 TUnfoldDensity.cxx:525
 TUnfoldDensity.cxx:526
 TUnfoldDensity.cxx:527
 TUnfoldDensity.cxx:528
 TUnfoldDensity.cxx:529
 TUnfoldDensity.cxx:530
 TUnfoldDensity.cxx:531
 TUnfoldDensity.cxx:532
 TUnfoldDensity.cxx:533
 TUnfoldDensity.cxx:534
 TUnfoldDensity.cxx:535
 TUnfoldDensity.cxx:536
 TUnfoldDensity.cxx:537
 TUnfoldDensity.cxx:538
 TUnfoldDensity.cxx:539
 TUnfoldDensity.cxx:540
 TUnfoldDensity.cxx:541
 TUnfoldDensity.cxx:542
 TUnfoldDensity.cxx:543
 TUnfoldDensity.cxx:544
 TUnfoldDensity.cxx:545
 TUnfoldDensity.cxx:546
 TUnfoldDensity.cxx:547
 TUnfoldDensity.cxx:548
 TUnfoldDensity.cxx:549
 TUnfoldDensity.cxx:550
 TUnfoldDensity.cxx:551
 TUnfoldDensity.cxx:552
 TUnfoldDensity.cxx:553
 TUnfoldDensity.cxx:554
 TUnfoldDensity.cxx:555
 TUnfoldDensity.cxx:556
 TUnfoldDensity.cxx:557
 TUnfoldDensity.cxx:558
 TUnfoldDensity.cxx:559
 TUnfoldDensity.cxx:560
 TUnfoldDensity.cxx:561
 TUnfoldDensity.cxx:562
 TUnfoldDensity.cxx:563
 TUnfoldDensity.cxx:564
 TUnfoldDensity.cxx:565
 TUnfoldDensity.cxx:566
 TUnfoldDensity.cxx:567
 TUnfoldDensity.cxx:568
 TUnfoldDensity.cxx:569
 TUnfoldDensity.cxx:570
 TUnfoldDensity.cxx:571
 TUnfoldDensity.cxx:572
 TUnfoldDensity.cxx:573
 TUnfoldDensity.cxx:574
 TUnfoldDensity.cxx:575
 TUnfoldDensity.cxx:576
 TUnfoldDensity.cxx:577
 TUnfoldDensity.cxx:578
 TUnfoldDensity.cxx:579
 TUnfoldDensity.cxx:580
 TUnfoldDensity.cxx:581
 TUnfoldDensity.cxx:582
 TUnfoldDensity.cxx:583
 TUnfoldDensity.cxx:584
 TUnfoldDensity.cxx:585
 TUnfoldDensity.cxx:586
 TUnfoldDensity.cxx:587
 TUnfoldDensity.cxx:588
 TUnfoldDensity.cxx:589
 TUnfoldDensity.cxx:590
 TUnfoldDensity.cxx:591
 TUnfoldDensity.cxx:592
 TUnfoldDensity.cxx:593
 TUnfoldDensity.cxx:594
 TUnfoldDensity.cxx:595
 TUnfoldDensity.cxx:596
 TUnfoldDensity.cxx:597
 TUnfoldDensity.cxx:598
 TUnfoldDensity.cxx:599
 TUnfoldDensity.cxx:600
 TUnfoldDensity.cxx:601
 TUnfoldDensity.cxx:602
 TUnfoldDensity.cxx:603
 TUnfoldDensity.cxx:604
 TUnfoldDensity.cxx:605
 TUnfoldDensity.cxx:606
 TUnfoldDensity.cxx:607
 TUnfoldDensity.cxx:608
 TUnfoldDensity.cxx:609
 TUnfoldDensity.cxx:610
 TUnfoldDensity.cxx:611
 TUnfoldDensity.cxx:612
 TUnfoldDensity.cxx:613
 TUnfoldDensity.cxx:614
 TUnfoldDensity.cxx:615
 TUnfoldDensity.cxx:616
 TUnfoldDensity.cxx:617
 TUnfoldDensity.cxx:618
 TUnfoldDensity.cxx:619
 TUnfoldDensity.cxx:620
 TUnfoldDensity.cxx:621
 TUnfoldDensity.cxx:622
 TUnfoldDensity.cxx:623
 TUnfoldDensity.cxx:624
 TUnfoldDensity.cxx:625
 TUnfoldDensity.cxx:626
 TUnfoldDensity.cxx:627
 TUnfoldDensity.cxx:628
 TUnfoldDensity.cxx:629
 TUnfoldDensity.cxx:630
 TUnfoldDensity.cxx:631
 TUnfoldDensity.cxx:632
 TUnfoldDensity.cxx:633
 TUnfoldDensity.cxx:634
 TUnfoldDensity.cxx:635
 TUnfoldDensity.cxx:636
 TUnfoldDensity.cxx:637
 TUnfoldDensity.cxx:638
 TUnfoldDensity.cxx:639
 TUnfoldDensity.cxx:640
 TUnfoldDensity.cxx:641
 TUnfoldDensity.cxx:642
 TUnfoldDensity.cxx:643
 TUnfoldDensity.cxx:644
 TUnfoldDensity.cxx:645
 TUnfoldDensity.cxx:646
 TUnfoldDensity.cxx:647
 TUnfoldDensity.cxx:648
 TUnfoldDensity.cxx:649
 TUnfoldDensity.cxx:650
 TUnfoldDensity.cxx:651
 TUnfoldDensity.cxx:652
 TUnfoldDensity.cxx:653
 TUnfoldDensity.cxx:654
 TUnfoldDensity.cxx:655
 TUnfoldDensity.cxx:656
 TUnfoldDensity.cxx:657
 TUnfoldDensity.cxx:658
 TUnfoldDensity.cxx:659
 TUnfoldDensity.cxx:660
 TUnfoldDensity.cxx:661
 TUnfoldDensity.cxx:662
 TUnfoldDensity.cxx:663
 TUnfoldDensity.cxx:664
 TUnfoldDensity.cxx:665
 TUnfoldDensity.cxx:666
 TUnfoldDensity.cxx:667
 TUnfoldDensity.cxx:668
 TUnfoldDensity.cxx:669
 TUnfoldDensity.cxx:670
 TUnfoldDensity.cxx:671
 TUnfoldDensity.cxx:672
 TUnfoldDensity.cxx:673
 TUnfoldDensity.cxx:674
 TUnfoldDensity.cxx:675
 TUnfoldDensity.cxx:676
 TUnfoldDensity.cxx:677
 TUnfoldDensity.cxx:678
 TUnfoldDensity.cxx:679
 TUnfoldDensity.cxx:680
 TUnfoldDensity.cxx:681
 TUnfoldDensity.cxx:682
 TUnfoldDensity.cxx:683
 TUnfoldDensity.cxx:684
 TUnfoldDensity.cxx:685
 TUnfoldDensity.cxx:686
 TUnfoldDensity.cxx:687
 TUnfoldDensity.cxx:688
 TUnfoldDensity.cxx:689
 TUnfoldDensity.cxx:690
 TUnfoldDensity.cxx:691
 TUnfoldDensity.cxx:692
 TUnfoldDensity.cxx:693
 TUnfoldDensity.cxx:694
 TUnfoldDensity.cxx:695
 TUnfoldDensity.cxx:696
 TUnfoldDensity.cxx:697
 TUnfoldDensity.cxx:698
 TUnfoldDensity.cxx:699
 TUnfoldDensity.cxx:700
 TUnfoldDensity.cxx:701
 TUnfoldDensity.cxx:702
 TUnfoldDensity.cxx:703
 TUnfoldDensity.cxx:704
 TUnfoldDensity.cxx:705
 TUnfoldDensity.cxx:706
 TUnfoldDensity.cxx:707
 TUnfoldDensity.cxx:708
 TUnfoldDensity.cxx:709
 TUnfoldDensity.cxx:710
 TUnfoldDensity.cxx:711
 TUnfoldDensity.cxx:712
 TUnfoldDensity.cxx:713
 TUnfoldDensity.cxx:714
 TUnfoldDensity.cxx:715
 TUnfoldDensity.cxx:716
 TUnfoldDensity.cxx:717
 TUnfoldDensity.cxx:718
 TUnfoldDensity.cxx:719
 TUnfoldDensity.cxx:720
 TUnfoldDensity.cxx:721
 TUnfoldDensity.cxx:722
 TUnfoldDensity.cxx:723
 TUnfoldDensity.cxx:724
 TUnfoldDensity.cxx:725
 TUnfoldDensity.cxx:726
 TUnfoldDensity.cxx:727
 TUnfoldDensity.cxx:728
 TUnfoldDensity.cxx:729
 TUnfoldDensity.cxx:730
 TUnfoldDensity.cxx:731
 TUnfoldDensity.cxx:732
 TUnfoldDensity.cxx:733
 TUnfoldDensity.cxx:734
 TUnfoldDensity.cxx:735
 TUnfoldDensity.cxx:736
 TUnfoldDensity.cxx:737
 TUnfoldDensity.cxx:738
 TUnfoldDensity.cxx:739
 TUnfoldDensity.cxx:740
 TUnfoldDensity.cxx:741
 TUnfoldDensity.cxx:742
 TUnfoldDensity.cxx:743
 TUnfoldDensity.cxx:744
 TUnfoldDensity.cxx:745
 TUnfoldDensity.cxx:746
 TUnfoldDensity.cxx:747
 TUnfoldDensity.cxx:748
 TUnfoldDensity.cxx:749
 TUnfoldDensity.cxx:750
 TUnfoldDensity.cxx:751
 TUnfoldDensity.cxx:752
 TUnfoldDensity.cxx:753
 TUnfoldDensity.cxx:754
 TUnfoldDensity.cxx:755
 TUnfoldDensity.cxx:756
 TUnfoldDensity.cxx:757
 TUnfoldDensity.cxx:758
 TUnfoldDensity.cxx:759
 TUnfoldDensity.cxx:760
 TUnfoldDensity.cxx:761
 TUnfoldDensity.cxx:762
 TUnfoldDensity.cxx:763
 TUnfoldDensity.cxx:764
 TUnfoldDensity.cxx:765
 TUnfoldDensity.cxx:766
 TUnfoldDensity.cxx:767
 TUnfoldDensity.cxx:768
 TUnfoldDensity.cxx:769
 TUnfoldDensity.cxx:770
 TUnfoldDensity.cxx:771
 TUnfoldDensity.cxx:772
 TUnfoldDensity.cxx:773
 TUnfoldDensity.cxx:774
 TUnfoldDensity.cxx:775
 TUnfoldDensity.cxx:776
 TUnfoldDensity.cxx:777
 TUnfoldDensity.cxx:778
 TUnfoldDensity.cxx:779
 TUnfoldDensity.cxx:780
 TUnfoldDensity.cxx:781
 TUnfoldDensity.cxx:782
 TUnfoldDensity.cxx:783
 TUnfoldDensity.cxx:784
 TUnfoldDensity.cxx:785
 TUnfoldDensity.cxx:786
 TUnfoldDensity.cxx:787
 TUnfoldDensity.cxx:788
 TUnfoldDensity.cxx:789
 TUnfoldDensity.cxx:790
 TUnfoldDensity.cxx:791
 TUnfoldDensity.cxx:792
 TUnfoldDensity.cxx:793
 TUnfoldDensity.cxx:794
 TUnfoldDensity.cxx:795
 TUnfoldDensity.cxx:796
 TUnfoldDensity.cxx:797
 TUnfoldDensity.cxx:798
 TUnfoldDensity.cxx:799
 TUnfoldDensity.cxx:800
 TUnfoldDensity.cxx:801
 TUnfoldDensity.cxx:802
 TUnfoldDensity.cxx:803
 TUnfoldDensity.cxx:804
 TUnfoldDensity.cxx:805
 TUnfoldDensity.cxx:806
 TUnfoldDensity.cxx:807
 TUnfoldDensity.cxx:808
 TUnfoldDensity.cxx:809
 TUnfoldDensity.cxx:810
 TUnfoldDensity.cxx:811
 TUnfoldDensity.cxx:812
 TUnfoldDensity.cxx:813
 TUnfoldDensity.cxx:814
 TUnfoldDensity.cxx:815
 TUnfoldDensity.cxx:816
 TUnfoldDensity.cxx:817
 TUnfoldDensity.cxx:818
 TUnfoldDensity.cxx:819
 TUnfoldDensity.cxx:820
 TUnfoldDensity.cxx:821
 TUnfoldDensity.cxx:822
 TUnfoldDensity.cxx:823
 TUnfoldDensity.cxx:824
 TUnfoldDensity.cxx:825
 TUnfoldDensity.cxx:826
 TUnfoldDensity.cxx:827
 TUnfoldDensity.cxx:828
 TUnfoldDensity.cxx:829
 TUnfoldDensity.cxx:830
 TUnfoldDensity.cxx:831
 TUnfoldDensity.cxx:832
 TUnfoldDensity.cxx:833
 TUnfoldDensity.cxx:834
 TUnfoldDensity.cxx:835
 TUnfoldDensity.cxx:836
 TUnfoldDensity.cxx:837
 TUnfoldDensity.cxx:838
 TUnfoldDensity.cxx:839
 TUnfoldDensity.cxx:840
 TUnfoldDensity.cxx:841
 TUnfoldDensity.cxx:842
 TUnfoldDensity.cxx:843
 TUnfoldDensity.cxx:844
 TUnfoldDensity.cxx:845
 TUnfoldDensity.cxx:846
 TUnfoldDensity.cxx:847
 TUnfoldDensity.cxx:848
 TUnfoldDensity.cxx:849
 TUnfoldDensity.cxx:850
 TUnfoldDensity.cxx:851
 TUnfoldDensity.cxx:852
 TUnfoldDensity.cxx:853
 TUnfoldDensity.cxx:854
 TUnfoldDensity.cxx:855
 TUnfoldDensity.cxx:856
 TUnfoldDensity.cxx:857
 TUnfoldDensity.cxx:858
 TUnfoldDensity.cxx:859
 TUnfoldDensity.cxx:860
 TUnfoldDensity.cxx:861
 TUnfoldDensity.cxx:862
 TUnfoldDensity.cxx:863
 TUnfoldDensity.cxx:864
 TUnfoldDensity.cxx:865
 TUnfoldDensity.cxx:866
 TUnfoldDensity.cxx:867
 TUnfoldDensity.cxx:868
 TUnfoldDensity.cxx:869
 TUnfoldDensity.cxx:870
 TUnfoldDensity.cxx:871
 TUnfoldDensity.cxx:872
 TUnfoldDensity.cxx:873
 TUnfoldDensity.cxx:874
 TUnfoldDensity.cxx:875
 TUnfoldDensity.cxx:876
 TUnfoldDensity.cxx:877
 TUnfoldDensity.cxx:878
 TUnfoldDensity.cxx:879
 TUnfoldDensity.cxx:880
 TUnfoldDensity.cxx:881
 TUnfoldDensity.cxx:882
 TUnfoldDensity.cxx:883
 TUnfoldDensity.cxx:884
 TUnfoldDensity.cxx:885
 TUnfoldDensity.cxx:886
 TUnfoldDensity.cxx:887
 TUnfoldDensity.cxx:888
 TUnfoldDensity.cxx:889
 TUnfoldDensity.cxx:890
 TUnfoldDensity.cxx:891
 TUnfoldDensity.cxx:892
 TUnfoldDensity.cxx:893
 TUnfoldDensity.cxx:894
 TUnfoldDensity.cxx:895
 TUnfoldDensity.cxx:896
 TUnfoldDensity.cxx:897
 TUnfoldDensity.cxx:898
 TUnfoldDensity.cxx:899
 TUnfoldDensity.cxx:900
 TUnfoldDensity.cxx:901
 TUnfoldDensity.cxx:902
 TUnfoldDensity.cxx:903
 TUnfoldDensity.cxx:904
 TUnfoldDensity.cxx:905
 TUnfoldDensity.cxx:906
 TUnfoldDensity.cxx:907
 TUnfoldDensity.cxx:908
 TUnfoldDensity.cxx:909
 TUnfoldDensity.cxx:910
 TUnfoldDensity.cxx:911
 TUnfoldDensity.cxx:912
 TUnfoldDensity.cxx:913
 TUnfoldDensity.cxx:914
 TUnfoldDensity.cxx:915
 TUnfoldDensity.cxx:916
 TUnfoldDensity.cxx:917
 TUnfoldDensity.cxx:918
 TUnfoldDensity.cxx:919
 TUnfoldDensity.cxx:920
 TUnfoldDensity.cxx:921
 TUnfoldDensity.cxx:922
 TUnfoldDensity.cxx:923
 TUnfoldDensity.cxx:924
 TUnfoldDensity.cxx:925
 TUnfoldDensity.cxx:926
 TUnfoldDensity.cxx:927
 TUnfoldDensity.cxx:928
 TUnfoldDensity.cxx:929
 TUnfoldDensity.cxx:930
 TUnfoldDensity.cxx:931
 TUnfoldDensity.cxx:932
 TUnfoldDensity.cxx:933
 TUnfoldDensity.cxx:934
 TUnfoldDensity.cxx:935
 TUnfoldDensity.cxx:936
 TUnfoldDensity.cxx:937
 TUnfoldDensity.cxx:938
 TUnfoldDensity.cxx:939
 TUnfoldDensity.cxx:940
 TUnfoldDensity.cxx:941
 TUnfoldDensity.cxx:942
 TUnfoldDensity.cxx:943
 TUnfoldDensity.cxx:944
 TUnfoldDensity.cxx:945
 TUnfoldDensity.cxx:946
 TUnfoldDensity.cxx:947
 TUnfoldDensity.cxx:948
 TUnfoldDensity.cxx:949
 TUnfoldDensity.cxx:950
 TUnfoldDensity.cxx:951
 TUnfoldDensity.cxx:952
 TUnfoldDensity.cxx:953
 TUnfoldDensity.cxx:954
 TUnfoldDensity.cxx:955
 TUnfoldDensity.cxx:956
 TUnfoldDensity.cxx:957
 TUnfoldDensity.cxx:958
 TUnfoldDensity.cxx:959
 TUnfoldDensity.cxx:960
 TUnfoldDensity.cxx:961
 TUnfoldDensity.cxx:962
 TUnfoldDensity.cxx:963
 TUnfoldDensity.cxx:964
 TUnfoldDensity.cxx:965
 TUnfoldDensity.cxx:966
 TUnfoldDensity.cxx:967
 TUnfoldDensity.cxx:968
 TUnfoldDensity.cxx:969
 TUnfoldDensity.cxx:970
 TUnfoldDensity.cxx:971
 TUnfoldDensity.cxx:972
 TUnfoldDensity.cxx:973
 TUnfoldDensity.cxx:974
 TUnfoldDensity.cxx:975
 TUnfoldDensity.cxx:976
 TUnfoldDensity.cxx:977
 TUnfoldDensity.cxx:978
 TUnfoldDensity.cxx:979
 TUnfoldDensity.cxx:980
 TUnfoldDensity.cxx:981
 TUnfoldDensity.cxx:982
 TUnfoldDensity.cxx:983
 TUnfoldDensity.cxx:984
 TUnfoldDensity.cxx:985
 TUnfoldDensity.cxx:986
 TUnfoldDensity.cxx:987
 TUnfoldDensity.cxx:988
 TUnfoldDensity.cxx:989
 TUnfoldDensity.cxx:990
 TUnfoldDensity.cxx:991
 TUnfoldDensity.cxx:992
 TUnfoldDensity.cxx:993
 TUnfoldDensity.cxx:994
 TUnfoldDensity.cxx:995
 TUnfoldDensity.cxx:996
 TUnfoldDensity.cxx:997
 TUnfoldDensity.cxx:998
 TUnfoldDensity.cxx:999
 TUnfoldDensity.cxx:1000
 TUnfoldDensity.cxx:1001
 TUnfoldDensity.cxx:1002
 TUnfoldDensity.cxx:1003
 TUnfoldDensity.cxx:1004
 TUnfoldDensity.cxx:1005
 TUnfoldDensity.cxx:1006
 TUnfoldDensity.cxx:1007
 TUnfoldDensity.cxx:1008
 TUnfoldDensity.cxx:1009
 TUnfoldDensity.cxx:1010
 TUnfoldDensity.cxx:1011
 TUnfoldDensity.cxx:1012
 TUnfoldDensity.cxx:1013
 TUnfoldDensity.cxx:1014
 TUnfoldDensity.cxx:1015
 TUnfoldDensity.cxx:1016
 TUnfoldDensity.cxx:1017
 TUnfoldDensity.cxx:1018
 TUnfoldDensity.cxx:1019
 TUnfoldDensity.cxx:1020
 TUnfoldDensity.cxx:1021
 TUnfoldDensity.cxx:1022
 TUnfoldDensity.cxx:1023
 TUnfoldDensity.cxx:1024
 TUnfoldDensity.cxx:1025
 TUnfoldDensity.cxx:1026
 TUnfoldDensity.cxx:1027
 TUnfoldDensity.cxx:1028
 TUnfoldDensity.cxx:1029
 TUnfoldDensity.cxx:1030
 TUnfoldDensity.cxx:1031
 TUnfoldDensity.cxx:1032
 TUnfoldDensity.cxx:1033
 TUnfoldDensity.cxx:1034
 TUnfoldDensity.cxx:1035
 TUnfoldDensity.cxx:1036
 TUnfoldDensity.cxx:1037
 TUnfoldDensity.cxx:1038
 TUnfoldDensity.cxx:1039
 TUnfoldDensity.cxx:1040
 TUnfoldDensity.cxx:1041
 TUnfoldDensity.cxx:1042
 TUnfoldDensity.cxx:1043
 TUnfoldDensity.cxx:1044
 TUnfoldDensity.cxx:1045
 TUnfoldDensity.cxx:1046
 TUnfoldDensity.cxx:1047
 TUnfoldDensity.cxx:1048
 TUnfoldDensity.cxx:1049
 TUnfoldDensity.cxx:1050
 TUnfoldDensity.cxx:1051
 TUnfoldDensity.cxx:1052
 TUnfoldDensity.cxx:1053
 TUnfoldDensity.cxx:1054
 TUnfoldDensity.cxx:1055
 TUnfoldDensity.cxx:1056
 TUnfoldDensity.cxx:1057
 TUnfoldDensity.cxx:1058
 TUnfoldDensity.cxx:1059
 TUnfoldDensity.cxx:1060
 TUnfoldDensity.cxx:1061
 TUnfoldDensity.cxx:1062
 TUnfoldDensity.cxx:1063
 TUnfoldDensity.cxx:1064
 TUnfoldDensity.cxx:1065
 TUnfoldDensity.cxx:1066
 TUnfoldDensity.cxx:1067
 TUnfoldDensity.cxx:1068
 TUnfoldDensity.cxx:1069
 TUnfoldDensity.cxx:1070
 TUnfoldDensity.cxx:1071
 TUnfoldDensity.cxx:1072
 TUnfoldDensity.cxx:1073
 TUnfoldDensity.cxx:1074
 TUnfoldDensity.cxx:1075
 TUnfoldDensity.cxx:1076
 TUnfoldDensity.cxx:1077
 TUnfoldDensity.cxx:1078
 TUnfoldDensity.cxx:1079
 TUnfoldDensity.cxx:1080
 TUnfoldDensity.cxx:1081
 TUnfoldDensity.cxx:1082
 TUnfoldDensity.cxx:1083
 TUnfoldDensity.cxx:1084
 TUnfoldDensity.cxx:1085
 TUnfoldDensity.cxx:1086
 TUnfoldDensity.cxx:1087
 TUnfoldDensity.cxx:1088
 TUnfoldDensity.cxx:1089
 TUnfoldDensity.cxx:1090
 TUnfoldDensity.cxx:1091
 TUnfoldDensity.cxx:1092
 TUnfoldDensity.cxx:1093
 TUnfoldDensity.cxx:1094
 TUnfoldDensity.cxx:1095
 TUnfoldDensity.cxx:1096
 TUnfoldDensity.cxx:1097
 TUnfoldDensity.cxx:1098
 TUnfoldDensity.cxx:1099
 TUnfoldDensity.cxx:1100
 TUnfoldDensity.cxx:1101
 TUnfoldDensity.cxx:1102
 TUnfoldDensity.cxx:1103
 TUnfoldDensity.cxx:1104
 TUnfoldDensity.cxx:1105
 TUnfoldDensity.cxx:1106
 TUnfoldDensity.cxx:1107
 TUnfoldDensity.cxx:1108
 TUnfoldDensity.cxx:1109
 TUnfoldDensity.cxx:1110
 TUnfoldDensity.cxx:1111
 TUnfoldDensity.cxx:1112
 TUnfoldDensity.cxx:1113
 TUnfoldDensity.cxx:1114
 TUnfoldDensity.cxx:1115
 TUnfoldDensity.cxx:1116
 TUnfoldDensity.cxx:1117
 TUnfoldDensity.cxx:1118
 TUnfoldDensity.cxx:1119
 TUnfoldDensity.cxx:1120
 TUnfoldDensity.cxx:1121
 TUnfoldDensity.cxx:1122
 TUnfoldDensity.cxx:1123
 TUnfoldDensity.cxx:1124
 TUnfoldDensity.cxx:1125
 TUnfoldDensity.cxx:1126
 TUnfoldDensity.cxx:1127
 TUnfoldDensity.cxx:1128
 TUnfoldDensity.cxx:1129
 TUnfoldDensity.cxx:1130
 TUnfoldDensity.cxx:1131
 TUnfoldDensity.cxx:1132
 TUnfoldDensity.cxx:1133
 TUnfoldDensity.cxx:1134
 TUnfoldDensity.cxx:1135
 TUnfoldDensity.cxx:1136
 TUnfoldDensity.cxx:1137
 TUnfoldDensity.cxx:1138
 TUnfoldDensity.cxx:1139
 TUnfoldDensity.cxx:1140
 TUnfoldDensity.cxx:1141
 TUnfoldDensity.cxx:1142
 TUnfoldDensity.cxx:1143
 TUnfoldDensity.cxx:1144
 TUnfoldDensity.cxx:1145
 TUnfoldDensity.cxx:1146
 TUnfoldDensity.cxx:1147
 TUnfoldDensity.cxx:1148
 TUnfoldDensity.cxx:1149
 TUnfoldDensity.cxx:1150
 TUnfoldDensity.cxx:1151
 TUnfoldDensity.cxx:1152
 TUnfoldDensity.cxx:1153
 TUnfoldDensity.cxx:1154
 TUnfoldDensity.cxx:1155
 TUnfoldDensity.cxx:1156
 TUnfoldDensity.cxx:1157
 TUnfoldDensity.cxx:1158
 TUnfoldDensity.cxx:1159
 TUnfoldDensity.cxx:1160
 TUnfoldDensity.cxx:1161
 TUnfoldDensity.cxx:1162
 TUnfoldDensity.cxx:1163
 TUnfoldDensity.cxx:1164
 TUnfoldDensity.cxx:1165
 TUnfoldDensity.cxx:1166
 TUnfoldDensity.cxx:1167
 TUnfoldDensity.cxx:1168
 TUnfoldDensity.cxx:1169
 TUnfoldDensity.cxx:1170
 TUnfoldDensity.cxx:1171
 TUnfoldDensity.cxx:1172
 TUnfoldDensity.cxx:1173
 TUnfoldDensity.cxx:1174
 TUnfoldDensity.cxx:1175
 TUnfoldDensity.cxx:1176
 TUnfoldDensity.cxx:1177
 TUnfoldDensity.cxx:1178
 TUnfoldDensity.cxx:1179
 TUnfoldDensity.cxx:1180
 TUnfoldDensity.cxx:1181
 TUnfoldDensity.cxx:1182
 TUnfoldDensity.cxx:1183
 TUnfoldDensity.cxx:1184
 TUnfoldDensity.cxx:1185
 TUnfoldDensity.cxx:1186
 TUnfoldDensity.cxx:1187
 TUnfoldDensity.cxx:1188
 TUnfoldDensity.cxx:1189
 TUnfoldDensity.cxx:1190
 TUnfoldDensity.cxx:1191
 TUnfoldDensity.cxx:1192
 TUnfoldDensity.cxx:1193
 TUnfoldDensity.cxx:1194
 TUnfoldDensity.cxx:1195
 TUnfoldDensity.cxx:1196
 TUnfoldDensity.cxx:1197
 TUnfoldDensity.cxx:1198
 TUnfoldDensity.cxx:1199
 TUnfoldDensity.cxx:1200
 TUnfoldDensity.cxx:1201
 TUnfoldDensity.cxx:1202
 TUnfoldDensity.cxx:1203
 TUnfoldDensity.cxx:1204
 TUnfoldDensity.cxx:1205
 TUnfoldDensity.cxx:1206
 TUnfoldDensity.cxx:1207
 TUnfoldDensity.cxx:1208
 TUnfoldDensity.cxx:1209
 TUnfoldDensity.cxx:1210
 TUnfoldDensity.cxx:1211
 TUnfoldDensity.cxx:1212
 TUnfoldDensity.cxx:1213
 TUnfoldDensity.cxx:1214
 TUnfoldDensity.cxx:1215
 TUnfoldDensity.cxx:1216
 TUnfoldDensity.cxx:1217
 TUnfoldDensity.cxx:1218
 TUnfoldDensity.cxx:1219
 TUnfoldDensity.cxx:1220
 TUnfoldDensity.cxx:1221
 TUnfoldDensity.cxx:1222
 TUnfoldDensity.cxx:1223
 TUnfoldDensity.cxx:1224
 TUnfoldDensity.cxx:1225
 TUnfoldDensity.cxx:1226
 TUnfoldDensity.cxx:1227
 TUnfoldDensity.cxx:1228
 TUnfoldDensity.cxx:1229
 TUnfoldDensity.cxx:1230
 TUnfoldDensity.cxx:1231
 TUnfoldDensity.cxx:1232
 TUnfoldDensity.cxx:1233
 TUnfoldDensity.cxx:1234
 TUnfoldDensity.cxx:1235
 TUnfoldDensity.cxx:1236
 TUnfoldDensity.cxx:1237
 TUnfoldDensity.cxx:1238
 TUnfoldDensity.cxx:1239
 TUnfoldDensity.cxx:1240
 TUnfoldDensity.cxx:1241
 TUnfoldDensity.cxx:1242
 TUnfoldDensity.cxx:1243
 TUnfoldDensity.cxx:1244
 TUnfoldDensity.cxx:1245
 TUnfoldDensity.cxx:1246
 TUnfoldDensity.cxx:1247
 TUnfoldDensity.cxx:1248
 TUnfoldDensity.cxx:1249
 TUnfoldDensity.cxx:1250
 TUnfoldDensity.cxx:1251
 TUnfoldDensity.cxx:1252
 TUnfoldDensity.cxx:1253
 TUnfoldDensity.cxx:1254
 TUnfoldDensity.cxx:1255
 TUnfoldDensity.cxx:1256
 TUnfoldDensity.cxx:1257
 TUnfoldDensity.cxx:1258
 TUnfoldDensity.cxx:1259
 TUnfoldDensity.cxx:1260
 TUnfoldDensity.cxx:1261
 TUnfoldDensity.cxx:1262
 TUnfoldDensity.cxx:1263
 TUnfoldDensity.cxx:1264
 TUnfoldDensity.cxx:1265
 TUnfoldDensity.cxx:1266
 TUnfoldDensity.cxx:1267
 TUnfoldDensity.cxx:1268
 TUnfoldDensity.cxx:1269
 TUnfoldDensity.cxx:1270
 TUnfoldDensity.cxx:1271
 TUnfoldDensity.cxx:1272
 TUnfoldDensity.cxx:1273
 TUnfoldDensity.cxx:1274
 TUnfoldDensity.cxx:1275
 TUnfoldDensity.cxx:1276
 TUnfoldDensity.cxx:1277
 TUnfoldDensity.cxx:1278
 TUnfoldDensity.cxx:1279
 TUnfoldDensity.cxx:1280
 TUnfoldDensity.cxx:1281
 TUnfoldDensity.cxx:1282
 TUnfoldDensity.cxx:1283
 TUnfoldDensity.cxx:1284
 TUnfoldDensity.cxx:1285
 TUnfoldDensity.cxx:1286
 TUnfoldDensity.cxx:1287
 TUnfoldDensity.cxx:1288
 TUnfoldDensity.cxx:1289
 TUnfoldDensity.cxx:1290
 TUnfoldDensity.cxx:1291
 TUnfoldDensity.cxx:1292
 TUnfoldDensity.cxx:1293
 TUnfoldDensity.cxx:1294
 TUnfoldDensity.cxx:1295
 TUnfoldDensity.cxx:1296
 TUnfoldDensity.cxx:1297
 TUnfoldDensity.cxx:1298
 TUnfoldDensity.cxx:1299
 TUnfoldDensity.cxx:1300
 TUnfoldDensity.cxx:1301
 TUnfoldDensity.cxx:1302
 TUnfoldDensity.cxx:1303
 TUnfoldDensity.cxx:1304
 TUnfoldDensity.cxx:1305
 TUnfoldDensity.cxx:1306
 TUnfoldDensity.cxx:1307
 TUnfoldDensity.cxx:1308
 TUnfoldDensity.cxx:1309
 TUnfoldDensity.cxx:1310
 TUnfoldDensity.cxx:1311
 TUnfoldDensity.cxx:1312
 TUnfoldDensity.cxx:1313
 TUnfoldDensity.cxx:1314
 TUnfoldDensity.cxx:1315
 TUnfoldDensity.cxx:1316
 TUnfoldDensity.cxx:1317
 TUnfoldDensity.cxx:1318
 TUnfoldDensity.cxx:1319
 TUnfoldDensity.cxx:1320
 TUnfoldDensity.cxx:1321
 TUnfoldDensity.cxx:1322
 TUnfoldDensity.cxx:1323
 TUnfoldDensity.cxx:1324
 TUnfoldDensity.cxx:1325
 TUnfoldDensity.cxx:1326
 TUnfoldDensity.cxx:1327
 TUnfoldDensity.cxx:1328
 TUnfoldDensity.cxx:1329
 TUnfoldDensity.cxx:1330
 TUnfoldDensity.cxx:1331
 TUnfoldDensity.cxx:1332
 TUnfoldDensity.cxx:1333
 TUnfoldDensity.cxx:1334
 TUnfoldDensity.cxx:1335
 TUnfoldDensity.cxx:1336
 TUnfoldDensity.cxx:1337
 TUnfoldDensity.cxx:1338
 TUnfoldDensity.cxx:1339
 TUnfoldDensity.cxx:1340
 TUnfoldDensity.cxx:1341
 TUnfoldDensity.cxx:1342
 TUnfoldDensity.cxx:1343
 TUnfoldDensity.cxx:1344
 TUnfoldDensity.cxx:1345
 TUnfoldDensity.cxx:1346
 TUnfoldDensity.cxx:1347
 TUnfoldDensity.cxx:1348
 TUnfoldDensity.cxx:1349
 TUnfoldDensity.cxx:1350
 TUnfoldDensity.cxx:1351
 TUnfoldDensity.cxx:1352
 TUnfoldDensity.cxx:1353
 TUnfoldDensity.cxx:1354
 TUnfoldDensity.cxx:1355
 TUnfoldDensity.cxx:1356
 TUnfoldDensity.cxx:1357
 TUnfoldDensity.cxx:1358
 TUnfoldDensity.cxx:1359
 TUnfoldDensity.cxx:1360
 TUnfoldDensity.cxx:1361
 TUnfoldDensity.cxx:1362
 TUnfoldDensity.cxx:1363
 TUnfoldDensity.cxx:1364
 TUnfoldDensity.cxx:1365
 TUnfoldDensity.cxx:1366
 TUnfoldDensity.cxx:1367
 TUnfoldDensity.cxx:1368
 TUnfoldDensity.cxx:1369
 TUnfoldDensity.cxx:1370
 TUnfoldDensity.cxx:1371
 TUnfoldDensity.cxx:1372
 TUnfoldDensity.cxx:1373
 TUnfoldDensity.cxx:1374
 TUnfoldDensity.cxx:1375
 TUnfoldDensity.cxx:1376
 TUnfoldDensity.cxx:1377
 TUnfoldDensity.cxx:1378
 TUnfoldDensity.cxx:1379
 TUnfoldDensity.cxx:1380
 TUnfoldDensity.cxx:1381
 TUnfoldDensity.cxx:1382
 TUnfoldDensity.cxx:1383
 TUnfoldDensity.cxx:1384
 TUnfoldDensity.cxx:1385
 TUnfoldDensity.cxx:1386
 TUnfoldDensity.cxx:1387
 TUnfoldDensity.cxx:1388
 TUnfoldDensity.cxx:1389
 TUnfoldDensity.cxx:1390
 TUnfoldDensity.cxx:1391
 TUnfoldDensity.cxx:1392
 TUnfoldDensity.cxx:1393
 TUnfoldDensity.cxx:1394
 TUnfoldDensity.cxx:1395
 TUnfoldDensity.cxx:1396
 TUnfoldDensity.cxx:1397
 TUnfoldDensity.cxx:1398
 TUnfoldDensity.cxx:1399
 TUnfoldDensity.cxx:1400
 TUnfoldDensity.cxx:1401
 TUnfoldDensity.cxx:1402
 TUnfoldDensity.cxx:1403
 TUnfoldDensity.cxx:1404
 TUnfoldDensity.cxx:1405
 TUnfoldDensity.cxx:1406
 TUnfoldDensity.cxx:1407
 TUnfoldDensity.cxx:1408
 TUnfoldDensity.cxx:1409
 TUnfoldDensity.cxx:1410
 TUnfoldDensity.cxx:1411
 TUnfoldDensity.cxx:1412
 TUnfoldDensity.cxx:1413
 TUnfoldDensity.cxx:1414
 TUnfoldDensity.cxx:1415
 TUnfoldDensity.cxx:1416
 TUnfoldDensity.cxx:1417
 TUnfoldDensity.cxx:1418
 TUnfoldDensity.cxx:1419
 TUnfoldDensity.cxx:1420
 TUnfoldDensity.cxx:1421
 TUnfoldDensity.cxx:1422
 TUnfoldDensity.cxx:1423
 TUnfoldDensity.cxx:1424
 TUnfoldDensity.cxx:1425
 TUnfoldDensity.cxx:1426
 TUnfoldDensity.cxx:1427
 TUnfoldDensity.cxx:1428
 TUnfoldDensity.cxx:1429
 TUnfoldDensity.cxx:1430
 TUnfoldDensity.cxx:1431
 TUnfoldDensity.cxx:1432
 TUnfoldDensity.cxx:1433
 TUnfoldDensity.cxx:1434
 TUnfoldDensity.cxx:1435
 TUnfoldDensity.cxx:1436
 TUnfoldDensity.cxx:1437
 TUnfoldDensity.cxx:1438
 TUnfoldDensity.cxx:1439
 TUnfoldDensity.cxx:1440
 TUnfoldDensity.cxx:1441
 TUnfoldDensity.cxx:1442
 TUnfoldDensity.cxx:1443
 TUnfoldDensity.cxx:1444
 TUnfoldDensity.cxx:1445
 TUnfoldDensity.cxx:1446
 TUnfoldDensity.cxx:1447
 TUnfoldDensity.cxx:1448
 TUnfoldDensity.cxx:1449
 TUnfoldDensity.cxx:1450
 TUnfoldDensity.cxx:1451
 TUnfoldDensity.cxx:1452
 TUnfoldDensity.cxx:1453
 TUnfoldDensity.cxx:1454
 TUnfoldDensity.cxx:1455
 TUnfoldDensity.cxx:1456
 TUnfoldDensity.cxx:1457
 TUnfoldDensity.cxx:1458
 TUnfoldDensity.cxx:1459
 TUnfoldDensity.cxx:1460
 TUnfoldDensity.cxx:1461
 TUnfoldDensity.cxx:1462
 TUnfoldDensity.cxx:1463
 TUnfoldDensity.cxx:1464
 TUnfoldDensity.cxx:1465
 TUnfoldDensity.cxx:1466
 TUnfoldDensity.cxx:1467
 TUnfoldDensity.cxx:1468
 TUnfoldDensity.cxx:1469
 TUnfoldDensity.cxx:1470
 TUnfoldDensity.cxx:1471
 TUnfoldDensity.cxx:1472
 TUnfoldDensity.cxx:1473
 TUnfoldDensity.cxx:1474
 TUnfoldDensity.cxx:1475
 TUnfoldDensity.cxx:1476
 TUnfoldDensity.cxx:1477
 TUnfoldDensity.cxx:1478
 TUnfoldDensity.cxx:1479
 TUnfoldDensity.cxx:1480
 TUnfoldDensity.cxx:1481
 TUnfoldDensity.cxx:1482
 TUnfoldDensity.cxx:1483
 TUnfoldDensity.cxx:1484
 TUnfoldDensity.cxx:1485
 TUnfoldDensity.cxx:1486
 TUnfoldDensity.cxx:1487
 TUnfoldDensity.cxx:1488
 TUnfoldDensity.cxx:1489
 TUnfoldDensity.cxx:1490
 TUnfoldDensity.cxx:1491
 TUnfoldDensity.cxx:1492
 TUnfoldDensity.cxx:1493
 TUnfoldDensity.cxx:1494
 TUnfoldDensity.cxx:1495
 TUnfoldDensity.cxx:1496
 TUnfoldDensity.cxx:1497
 TUnfoldDensity.cxx:1498
 TUnfoldDensity.cxx:1499
 TUnfoldDensity.cxx:1500
 TUnfoldDensity.cxx:1501
 TUnfoldDensity.cxx:1502
 TUnfoldDensity.cxx:1503
 TUnfoldDensity.cxx:1504