////////////////////////////////////////////////////////////////////////////
//                                                                        //
// This program should be run under root :                                //
//      root fluxunfold.C++                                               //
//                                                                        //
// Author(s) : T. Bretz  02/2002 <mailto:tbretz@astro.uni-wuerzburg.de>   //
//             W. Wittek 09/2002 <mailto:wittek@mppmu.mpg.de>             //
//                                                                        //
// this macro is prepared to be used in the analysis :                    //
//                                                                        //
//      the unfolding should be called by                                 //
//      doUnfolding(TH2D &tobeunfolded,      // (E-est, Theta)            //
//                  TH3D &migrationmatrix,   // (E-est, E-true, Theta)    //
//                  TH2D &unfolded)          // (E-true,Theta)            //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

#include <TMath.h>
#include <TRandom3.h>
#include <TVector.h>
#include <TMatrixD.h>
#include <TMatrix.h>
#include <TH1.h>
#include <TH2.h>
#include <TH3.h>
#include <TProfile.h>
#include <TF1.h>
#include <iostream.h>
#include <TMinuit.h>
#include <TCanvas.h>
#include <TMarker.h>

#include <fstream.h>
#include <iomanip.h>

TH1 *DrawMatrixClone(const TMatrixD &m, Option_t *opt="")
{
    const Int_t nrows = m.GetNrows();
    const Int_t ncols = m.GetNcols();

    TMatrix m2(nrows, ncols);
    for (int i=0; i<nrows; i++)
        for (int j=0; j<ncols; j++)
            m2(i, j) = m(i, j);

    TH2F *hist = new TH2F(m2);
    hist->SetBit(kCanDelete);
    hist->Draw(opt);
    hist->SetDirectory(NULL);

    return hist;

}

TH1 *DrawMatrixColClone(const TMatrixD &m, Option_t *opt="", Int_t col=0)
{
    const Int_t nrows = m.GetNrows();

    TVector vec(nrows);
    for (int i=0; i<nrows; i++)
        vec(i) = m(i, col);

    TH1F *hist = new TH1F("TVector","",nrows,0,nrows);
    for (int i=0; i<nrows; i++)
    {
      hist->SetBinContent(i+1, vec(i));
    }

    hist->SetBit(kCanDelete);
    hist->Draw(opt);
    hist->SetDirectory(NULL);

    return hist;
}


void PrintTH3Content(const TH3 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
    {
      for (Int_t j=1; j<=hist.GetNbinsY(); j++)
        for (Int_t k=1; k<=hist.GetNbinsZ(); k++)
            cout << hist.GetBinContent(i,j,k) << " \t";
      cout << endl << endl;
    }
}

void PrintTH3Error(const TH3 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
    {
      for (Int_t j=1; j<=hist.GetNbinsY(); j++)
        for (Int_t k=1; k<=hist.GetNbinsZ(); k++)
            cout << hist.GetBinError(i, j, k) << " \t";
      cout << endl << endl;
    }
}

void PrintTH2Content(const TH2 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
        for (Int_t j=1; j<=hist.GetNbinsY(); j++)
            cout << hist.GetBinContent(i,j) << " \t";
        cout << endl << endl;
}

void PrintTH2Error(const TH2 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
        for (Int_t j=1; j<=hist.GetNbinsY(); j++)
            cout << hist.GetBinError(i, j) << " \t";
        cout << endl << endl;
}

void PrintTH1Content(const TH1 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
        cout << hist.GetBinContent(i) << " \t";
    cout << endl << endl;
}

void PrintTH1Error(const TH1 &hist)
{
    cout << hist.GetName() << ": " << hist.GetTitle() << " <error>" << endl;
    cout << "-----------------------------------------------------" << endl;
    for (Int_t i=1; i<=hist.GetNbinsX(); i++)
        cout << hist.GetBinError(i) << " \t";
    cout << endl << endl;
}

void CopyCol(TMatrixD &m, const TH1 &h, Int_t col=0)
{
    const Int_t n = m.GetNrows();

    for (Int_t i=0; i<n; i++)
        m(i, col) = h.GetBinContent(i+1);
}

void CopyCol(TH1 &h, const TMatrixD &m, Int_t col=0)
{
    const Int_t n = m.GetNrows();

    for (Int_t i=0; i<n; i++)
        h.SetBinContent(i+1, m(i, col));
}

void CopyH2M(TMatrixD &m, const TH2 &h)
{
    const Int_t nx = m.GetNrows();
    const Int_t ny = m.GetNcols();

    for (Int_t i=0; i<nx; i++)
        for (Int_t j=0; j<ny; j++)
            m(i, j) = h.GetBinContent(i+1, j+1);
}

void CopySqr(TMatrixD &m, const TH1 &h)
{
    const Int_t nx = m.GetNrows();
    const Int_t ny = m.GetNcols();

    for (Int_t i=0; i<nx; i++)
        for (Int_t j=0; j<ny; j++)
        {
            const Double_t bin =  h.GetBinContent(i+1, j+1);
            m(i, j) = bin*bin;
        }
}

Double_t GetMatrixSumRow(const TMatrixD &m, Int_t row)
{
    const Int_t n = m.GetNcols();

    Double_t sum = 0;
    for (Int_t i=0; i<n; i++)
        sum += m(row, i);

    return sum;
}

Double_t GetMatrixSumDiag(const TMatrixD &m)
{
    const Int_t n = m.GetNcols();

    Double_t sum = 0;
    for (Int_t i=0; i<n; i++)
        sum += m(i, i);

    return sum;
}

Double_t GetMatrixSumCol(const TMatrixD &m, Int_t col=0)
{
    const Int_t n = m.GetNrows();

    Double_t sum = 0;
    for (Int_t i=0; i<n; i++)
        sum += m(i, col);

    return sum;
}
Double_t GetMatrixSum(const TMatrixD &m)
{
    const Int_t n = m.GetNrows();

    Double_t sum = 0;
    for (Int_t i=0; i<n; i++)
        sum += GetMatrixSumRow(m, i);

    return sum;
}

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// fcnSmooth   (used by SmoothMigrationMatrix)                                 //
//                                                                        //
// is called by MINUIT                                                    //
// for given values of the parameters it calculates                       //
//                     the function to be minimized                       //  
//                                                                        //
////////////////////////////////////////////////////////////////////////////
void fcnSmooth(Int_t &npar, Double_t *gin, Double_t &f, 
               Double_t *par, Int_t iflag);



////////////////////////////////////////////////////////////////////////////
//                                                                        //
// fcnTikhonov2   (used by Tikhonov2)                                     //
//                                                                        //
// is called by MINUIT                                                    //
// for given values of the parameters it calculates                       //
//                     the function to be minimized                       //  
//                                                                        //
////////////////////////////////////////////////////////////////////////////
void fcnTikhonov2(Int_t &npar, Double_t *gin, Double_t &f, 
                              Double_t *par, Int_t iflag);

////////////////////////////////////////////////////////////////////////////
//                                                                        //
// MUnfold                                                                //
//                                                                        //
// class for unfolding a 1-dimensional distribution                       //
//                                                                        //
// the methods used are described in :                                    //
//                                                                        //
//     V.B.Anykeyev et al., NIM A303 (1991) 350                           //
//     M. Schmelling, Nucl. Instr. and Meth. A 340 (1994) 400             //
//     M. Schmelling : "Numerische Methoden der Datenanalyse"             //
//                    Heidelberg, Maerz 1998                              //
//     M.Bertero, INFN/TC-88/2 (1988)                                     //
//                                                                        //
////////////////////////////////////////////////////////////////////////////
class MUnfold : public TObject
{
public:
    TString   bintitle;

    UInt_t    fNa;        // Number of bins in the distribution to be unfolded
    UInt_t    fNb;        // Number of bins in the unfolded distribution

    TMatrixD  fMigrat;    // migration matrix                  (fNa, fNb)
    TMatrixD  fMigraterr2;// error**2 of migration matrix      (fNa, fNb)

    TMatrixD  fMigOrig;    // original migration matrix         (fNa, fNb)
    TMatrixD  fMigOrigerr2;// error**2 oforiginal migr. matrix  (fNa, fNb)

    TMatrixD  fMigSmoo;    // smoothed migration matrix M       (fNa, fNb)
    TMatrixD  fMigSmooerr2;// error**2 of smoothed migr. matrix (fNa, fNb)
    TMatrixD  fMigChi2;    // chi2 contributions for smoothing  (fNa, fNb)

    TMatrixD  fVa;        // distribution to be unfolded       (fNa)
    TMatrixD  fVacov;     // error matrix of fVa               (fNa, fNa)
    TMatrixD  fVacovInv; // inverse of fVacov                 (fNa, fNa)
    Double_t  fSpurVacov; // Spur of fVacov

    //    UInt_t    fVaevents;  // total number of events
    UInt_t    fVapoints;  // number of significant measurements

    TMatrixD  fVb;        // unfolded distribution             (fNb)
    TMatrixD  fVbcov;     // error matrix of fVb               (fNb, fNb)

    TMatrixD  fVEps;      // prior distribution                (fNb)
    TMatrixDColumn fVEps0;

    Double_t  fW;         // weight
    Double_t  fWbest;     // best weight
    Int_t     ixbest;

    TMatrixD  fResult;    // unfolded distribution and errors  (fNb, 5)
    TMatrixD  fChi2;      // chisquared contribution           (fNa, 1)

    Double_t  fChisq;     // total chisquared
    Double_t  fNdf;       // number of degrees of freedom
    Double_t  fProb;      // chisquared probability

    TMatrixD  G;          // G = M * M(transposed)             (fNa, fNa)
    TVectorD  EigenValue; // vector of eigenvalues lambda of G (fNa)
    TMatrixD  Eigen;      // matrix of eigen vectors of G      (fNa, fNa)
    Double_t  RankG;      // rank of G
    Double_t  tau;        // 1 / lambda_max
    Double_t  EpsLambda;

    // quantities stored for each weight :
    TVectorD SpSig;       // Spur of covariance matrix of fVbcov
    TVectorD SpAR;        // effective rank of G^tilde
    TVectorD chisq;       // chi squared (measures agreement between
    // fVa and the folded fVb)
    TVectorD SecDer;      // regularization term = sum of (2nd der.)**2
    TVectorD ZerDer;      // regularization term = sum of (fVb)**2
    TVectorD Entrop;      // regularization term = reduced cross-entropy
    TVectorD DAR2;        //
    TVectorD Dsqbar;      //

    Double_t SpurAR;
    Double_t SpurSigma;
    Double_t SecDeriv;
    Double_t ZerDeriv;
    Double_t Entropy;
    Double_t DiffAR2;
    Double_t Chisq;
    Double_t D2bar;

    TMatrixD Chi2;

    //

    // plots versus weight
    Int_t    Nix;
    Double_t xmin;
    Double_t xmax;
    Double_t dlogx;

    TH1D *hBchisq;
    TH1D *hBSpAR;
    TH1D *hBDSpAR;
    TH1D *hBSpSig;
    TH1D *hBDSpSig;
    TH1D *hBSecDeriv;
    TH1D *hBDSecDeriv;
    TH1D *hBZerDeriv;
    TH1D *hBDZerDeriv;
    TH1D *hBEntropy;
    TH1D *hBDEntropy;
    TH1D *hBDAR2;
    TH1D *hBD2bar;

    //
    TH1D *hEigen;

    // plots for the best solution
    TH2D *fhmig;
    TH2D *shmig;
    TH2D *shmigChi2;

    TH1D *fhb0;

    TH1D *fha;

    TH1D *hprior;

    TH1D *hb;

    Double_t CalcSpurSigma(TMatrixD &T, Double_t norm=1)
    {
        Double_t spursigma = 0;

        for (UInt_t a=0; a<fNb; a++)
        {
            for (UInt_t b=0; b<fNb; b++)
            {
                fVbcov(a,b) = 0;

                for (UInt_t c=0; c<fNa; c++)
                    for (UInt_t d=0; d<fNa; d++)
                        fVbcov(a,b) += T(a,d)*fVacov(d,c)*T(b,c);

                fVbcov(a,b) *= norm*norm;
            }
            spursigma += fVbcov(a,a);
        }

        return spursigma;
    }

public:
    // -----------------------------------------------------------------------
    //
    // Constructor
    //              copy histograms into matrices
    //
    MUnfold(TH1D &ha, TH2D &hacov, TH2D &hmig)
        : fVEps(hmig.GetYaxis()->GetNbins(),1), fVEps0(fVEps, 0)
    {
        // ha      is the distribution to be unfolded
        // hacov   is the covariance matrix of ha
        // hmig    is the migration matrix;
        //         this matrix will be used in the unfolding
        //         unless SmoothMigrationMatrix(*hmigrat) is called;
        //         in the latter case hmigrat is smoothed
        //         and the smoothed matrix is used in the unfolding

        // Eigen values of the matrix G, which are smaller than EpsLambda
        // will be considered as being zero
        EpsLambda = 1.e-10;
        fW = 0.0;

        fNa  = hmig.GetXaxis()->GetNbins();
        const Double_t alow = hmig.GetXaxis()->GetXmin();
        const Double_t aup  = hmig.GetXaxis()->GetXmax();

        fNb  = hmig.GetYaxis()->GetNbins();
        const Double_t blow = hmig.GetYaxis()->GetXmin();
        const Double_t bup  = hmig.GetYaxis()->GetXmax();


        UInt_t Na = ha.GetNbinsX();
        if (fNa != Na)
        {
            cout << "MUnfold::MUnfold : dimensions do not match,  fNa = ";
            cout << fNa << ",   Na = " << Na << endl;
        }

        cout << "MUnfold::MUnfold :" << endl;
        cout << "==================" << endl;
        cout << "   fNa = " << fNa << ",   fNb = " << fNb << endl;

        // ------------------------

        fVa.ResizeTo(fNa, 1);
        CopyCol(fVa, ha, 0);

        cout << "   fVa = ";

        for (UInt_t i=0; i<fNa; i++)
            cout << fVa(i,0) << " \t";
        cout << endl;

        Double_t vaevents = GetMatrixSumCol(fVa, 0);
        cout << "   Total number of events in fVa = " << vaevents << endl;

        // ------------------------

        fChi2.ResizeTo(fNa,1);
        Chi2.ResizeTo(fNa,1);

        // ------------------------

        fVacov.ResizeTo(fNa, fNa);
        fSpurVacov = 0;

        CopyH2M(fVacov, hacov);

        fVapoints = 0;
        for (UInt_t i=0; i<fNa; i++)
            if (fVa(i,0)>0 && fVacov(i,i)<fVa(i,0)*fVa(i,0))
                fVapoints++;

        fSpurVacov = GetMatrixSumDiag(fVacov);

        //cout << "MUnfold::MUnfold :   fVacov = " << endl;
        //cout << "==============================" << endl;
        //fVacov.Print();

        cout << "   Number of significant points in fVa = ";
        cout << fVapoints << endl;

        cout << "   Spur of fVacov = ";
        cout << fSpurVacov << endl;

        // ------------------------

        fVacovInv.ResizeTo(fNa, fNa);
        fVacovInv = fVacov;
        fVacovInv.InvertPosDef();

        //cout << "MUnfold::MUnfold :   fVacovInv = " << endl;
        //cout << "==================================" << endl;
        //fVacovInv.Print();

        // ------------------------
        // fMigrat is the migration matrix to be used in the unfolding;
        // fMigrat may be overwritten by SmoothMigrationMatrix

        fMigrat.ResizeTo(fNa, fNb); // row, col

        CopyH2M(fMigrat, hmig);


        // ------------------------

        fMigraterr2.ResizeTo(fNa, fNb); // row, col
        CopySqr(fMigraterr2, hmig);

        // normaxlize

        for (UInt_t j=0; j<fNb; j++)
        {
            const Double_t sum = GetMatrixSumCol(fMigrat, j);

            if (sum==0)
                continue;

            TMatrixDColumn col1(fMigrat, j);
            col1 *= 1./sum;

            TMatrixDColumn col2(fMigraterr2, j);
            col2 *= 1./(sum*sum);
        }

        //cout << "MUnfold::MUnfold :   fMigrat = " << endl;
        //cout << "===============================" << endl;
        //fMigrat.Print();

        //cout << "MUnfold::MUnfold :   fMigraterr2 = " << endl;
        //cout << "===================================" << endl;
        //fMigraterr2.Print();

        // ------------------------
        G.ResizeTo(fNa, fNa);
        EigenValue.ResizeTo(fNa);
        Eigen.ResizeTo(fNa, fNa);

        fMigOrig.ResizeTo(fNa, fNb); 
        fMigOrigerr2.ResizeTo(fNa, fNb); 

        fMigSmoo.ResizeTo    (fNa, fNb); 
        fMigSmooerr2.ResizeTo(fNa, fNb); 
        fMigChi2.ResizeTo    (fNa, fNb); 

        // ------------------------

        fVEps0 = 1./fNb;

        //cout << "MUnfold::MUnfold :   Default prior distribution fVEps = " << endl;
        //cout << "========================================================" << endl;
        //fVEps.Print();

        // ------------------------

        fVb.ResizeTo(fNb,1);
        fVbcov.ResizeTo(fNb,fNb);

        // ----------------------------------------------------
        // number and range of weights to be scanned
        Nix  = 30;
        xmin = 1.e-5;
        xmax = 1.e5;
        dlogx = (log10(xmax)-log10(xmin)) / Nix;

        SpSig.ResizeTo (Nix);
        SpAR.ResizeTo  (Nix);
        chisq.ResizeTo (Nix);
        SecDer.ResizeTo(Nix);
        ZerDer.ResizeTo(Nix);
        Entrop.ResizeTo(Nix);
        DAR2.ResizeTo  (Nix);
        Dsqbar.ResizeTo(Nix);

        //------------------------------------
        // plots as a function of the iteration  number

        hBchisq = new TH1D("Bchisq", "chisq",
                           Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBSpAR  = new TH1D("BSpAR", "SpurAR",
                           Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDSpAR  = new TH1D("BDSpAR", "Delta(SpurAR)",
                            Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBSpSig = new TH1D("BSpSig", "SpurSigma/SpurC",
                           Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDSpSig = new TH1D("BDSpSig", "Delta(SpurSigma/SpurC)",
                            Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBSecDeriv = new TH1D("BSecDeriv", "Second Derivative squared",
                              Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDSecDeriv = new TH1D("BDSecDeriv", "Delta(Second Derivative squared)",
                               Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBZerDeriv = new TH1D("BZerDeriv", "Zero Derivative squared",
                              Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDZerDeriv = new TH1D("BDZerDeriv", "Delta(Zero Derivative squared)",
                               Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBEntropy = new TH1D("BEntrop", "Entropy",
                             Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDEntropy = new TH1D("BDEntrop", "Delta(Entropy)",
                              Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBDAR2 = new TH1D("BDAR2", "norm(AR-AR+)",
                          Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        hBD2bar = new TH1D("BD2bar", "(b_unfolded-b_ideal)**2",
                           Nix, log10(xmin)-dlogx/2.0, log10(xmax)-dlogx/2.0 );

        //-------------------------------------
        // original migration matrix
        fhmig = new TH2D("fMigrat", "Migration matrix",
                         fNa, alow, aup, fNb, blow, bup);
        fhmig->Sumw2();

        //-------------------------------------
        // smoothed migration matrix
        shmig = new TH2D("sMigrat", "Smoothed migration matrix",
                         fNa, alow, aup, fNb, blow, bup);
        shmig->Sumw2();

        //-------------------------------------
        // chi2 contributions for smoothing of migration matrix
        shmigChi2 = new TH2D("sMigratChi2", "Chi2 contr. for smoothing",
                             fNa, alow, aup, fNb, blow, bup);

        //-------------------------------------
        // eigen values of matrix G = M * M(transposed)
        hEigen = new TH1D("Eigen", "Eigen values of M*MT",
                          fNa, 0.5, fNa+0.5);

        //------------------------------------
        // Ideal distribution
        
        fhb0 = new TH1D("fhb0", "Ideal distribution", fNb, blow, bup);
        fhb0->Sumw2();
        

        //------------------------------------
        // Distribution to be unfolded
        fha = new TH1D("fha", "Distribution to be unfolded", fNa, alow, aup);
        fha->Sumw2();

        //------------------------------------
        // Prior distribution
        hprior = new TH1D("Prior", "Prior distribution", fNb, blow, bup);

        //------------------------------------
        // Unfolded distribution
        hb = new TH1D("DataSp", "Unfolded distribution", fNb, blow, bup);
        hb->Sumw2();

    }

    // -----------------------------------------------------------------------
    //
    // Define prior distribution to be a constant
    //
    void SetPriorConstant()
    {
        fVEps0 = 1./fNb;

        CopyCol(*hprior, fVEps);

        //cout << "SetPriorConstant : Prior distribution fVEps = " << endl;
        //cout << "==============================================" << endl;
        //fVEps.Print();
    }

    // -----------------------------------------------------------------------
    //
    // Take prior distribution from the histogram 'ha'
    // which may have a different binning than 'hprior'
    //
    Bool_t SetPriorRebin(TH1D &ha)
    {
        // ------------------------------------------------------------------
        //
        // fill the contents of histogram 'ha' into the histogram 'hrior';
        // the histograms need not have the same binning;
        // if the binnings are different, the bin contents of histogram 'ha'
        //    are distributed properly (linearly) over the bins of 'hprior'
        //

        const Int_t    na   = ha.GetNbinsX();
        const Double_t alow = ha.GetBinLowEdge(1);
        const Double_t aup  = ha.GetBinLowEdge(na+1);

        const Int_t    nb   = hprior->GetNbinsX();
        const Double_t blow = hprior->GetBinLowEdge(1);
        const Double_t bup  = hprior->GetBinLowEdge(nb+1);

        // check whether there is an overlap
        //       between the x ranges of the 2 histograms
        if (alow>bup || aup<blow)
        {
            cout << "Rebinning not possible because there is no overlap of the x ranges of the two histograms" << endl;
            return kFALSE;
        }

        // there is an overlap
        //********************
        Double_t sum = 0;
        for (Int_t j=1; j<=nb; j++)
        {
            const Double_t yl = hprior->GetBinLowEdge(j);
            const Double_t yh = hprior->GetBinLowEdge(j+1);

            // search bins of histogram ha which contribute
            // to bin j of histogram hb
            //----------------
            Int_t il=0;
            Int_t ih=0;
            for (Int_t i=2; i<=na+1; i++)
            {
                const Double_t xl = ha.GetBinLowEdge(i);
                if (xl>yl)
                {
                    il = i-1;

                    //.................................
                    ih = 0;
                    for (Int_t k=(il+1); k<=(na+1); k++)
                    {
                        const Double_t xh = ha.GetBinLowEdge(k);
                        if (xh >= yh)
                        {
                            ih = k-1;
                            break;
                        }
                    }
                    //.................................
                    if (ih == 0)
                        ih = na;
                    break;
                }
            }
            //----------------
            if (il == 0)
            {
                cout << "Something is wrong " << endl;
                cout << "          na, alow, aup = " << na << ",  " << alow
                    << ",  " << aup << endl;
                cout << "          nb, blow, bup = " << nb << ",  " << blow
                    << ",  " << bup << endl;
                return kFALSE;
            }

            Double_t content=0;
            // sum up the contribution to bin j
            for (Int_t i=il; i<=ih; i++)
            {
                const Double_t xl = ha.GetBinLowEdge(i);
                const Double_t xh = ha.GetBinLowEdge(i+1);
                const Double_t bina = xh-xl;

                if (xl<yl  &&  xh<yh)
                    content += ha.GetBinContent(i) * (xh-yl) / bina;
                else
                    if (xl<yl  &&  xh>=yh)
                        content += ha.GetBinContent(i) * (yh-yl) / bina;
                    else
                        if (xl>=yl  &&  xh<yh)
                            content += ha.GetBinContent(i);
                        else if (xl>=yl  &&  xh>=yh)
                            content += ha.GetBinContent(i) * (yh-xl) / bina;
            }
            hprior->SetBinContent(j, content);
            sum += content;
        }

        // normalize histogram hb
        if (sum==0)
        {
            cout << "histogram hb is empty; sum of weights in ha = ";
            cout << ha.GetSumOfWeights() << endl;
            return kFALSE;
        }

        for (Int_t j=1; j<=nb; j++)
        {
            const Double_t content = hprior->GetBinContent(j)/sum;
            hprior->SetBinContent(j, content);
            fVEps0(j-1) = content;
        }

        //cout << "SetPriorRebin : Prior distribution fVEps = " << endl;
        //cout << "===========================================" << endl;
        //fVEps.Print();

        return kTRUE;
    }


    // -----------------------------------------------------------------------
    //
    // Set prior distribution to a given distribution 'hpr'
    //
    Bool_t SetPriorInput(TH1D &hpr)
    {
        CopyCol(fVEps, hpr);

        const Double_t sum = GetMatrixSumCol(fVEps, 0);

        if (sum<=0)
        {
            cout << "MUnfold::SetPriorInput: invalid prior distribution" << endl;
            return kFALSE;
        }

        // normalize prior distribution
        fVEps0 *= 1./sum;

        CopyCol(*hprior, fVEps);

        //cout << "SetPriorInput : Prior distribution fVEps = " << endl;
        //cout << "===========================================" << endl;
        //fVEps.Print();

        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // Define prior distribution to be a power law
    // use input distribution 'hprior' only
    //           for defining the histogram parameters
    //
    Bool_t SetPriorPower(Double_t gamma)
    {
        // generate distribution according to a power law
        //                        dN/dE = E^{-gamma}
        //  or with y = lo10(E),  E = 10^y :
        //                        dN/dy = ln10 * 10^{y*(1-gamma)}
        TH1D hpower(*hprior);

        const UInt_t   nbin = hprior->GetNbinsX();
        const Double_t xmin = hprior->GetBinLowEdge(1);
        const Double_t xmax = hprior->GetBinLowEdge(nbin+1);

        cout << "nbin, xmin, xmax = " << nbin << ",  ";
        cout << xmin << ",  " << xmax << endl;

        TF1* fpow = new TF1("fpow", "pow(10.0, x*(1.0-[0]))", xmin,xmax);
        fpow->SetParName  (0,"gamma");
        fpow->SetParameter(0, gamma );

        hpower.FillRandom("fpow", 100000);

        // fill prior distribution
        CopyCol(fVEps, hpower);

        const Double_t sum = GetMatrixSumCol(fVEps, 0);
        if (sum <= 0)
        {
            cout << "MUnfold::SetPriorPower : invalid prior distribution"  << endl;
            return kFALSE;
        }

        // normalize prior distribution
        fVEps0 *= 1./sum;
        CopyCol(*hprior, fVEps);

        //cout << "SetPriorPower : Prior distribution fVEps = " << endl;
        //cout << "===========================================" << endl;
        //fVEps.Print();

        return kTRUE;
    }


    // -----------------------------------------------------------------------
    //
    // Set the initial weight
    //
    Bool_t SetInitialWeight(Double_t &weight)
    {
        if (weight == 0.0)
        {
            TMatrixD v1(fVa, TMatrixD::kTransposeMult, fVacovInv);
            TMatrixD v2(v1, TMatrixD::kMult, fVa);
            weight = 1./sqrt(v2(0,0));
        }

        cout << "MUnfold::SetInitialWeight : Initial Weight = "
            << weight << endl;

        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // Print the unfolded distribution
    //
    void PrintResults()
    {
        cout << bintitle << endl;
        cout << "PrintResults : Unfolded distribution fResult " << endl;
        cout << "=============================================" << endl;
        //cout << "val, eparab, eplus, eminus, gcc = "  << endl;

        for (UInt_t i=0; i<fNb; i++)
        {
	  //    cout << fResult(i, 0) << " \t";
          //  cout << fResult(i, 1) << " \t";
          //  cout << fResult(i, 2) << " \t";
          //  cout << fResult(i, 3) << " \t";
          //  cout << fResult(i, 4) <<  endl;
        }
        cout << "Chisquared, NDF, chi2 probability, ixbest = "
            << fChisq << ",  "
            << fNdf << ",  " << fProb << ",  " << ixbest << endl;

    }


    // -----------------------------------------------------------------------
    //
    // Schmelling  : unfolding by minimizing the function Z
    //               by Gauss-Newton iteration
    //
    //               the weights are scanned between
    //               1.e-5*fWinitial and 1.e5*fWinitial
    //
    Bool_t Schmelling(TH1D &hb0)
    {
    
        //======================================================================
        // copy ideal distribution
        for (UInt_t i=1; i<=fNb; i++)
        {
            fhb0->SetBinContent(i, hb0.GetBinContent(i));
            fhb0->SetBinError  (i, hb0.GetBinError(i));
        }
    
        //-----------------------------------------------------------------------
        // Initialization
        // ==============

        Int_t numGiteration;
        Int_t MaxGiteration = 1000;

        TMatrixD alpha;
        alpha.ResizeTo(fNa, 1);


        //-----------------------------------------------------------------------
        // Newton iteration
        // ================

        Double_t dga2;
        Double_t dga2old;
        Double_t EpsG = 1.e-12;

        TMatrixD wZdp_inv(fNa, fNa);
        TMatrixD d(fNb, 1);
        TMatrixD p(fNb, 1);

        TMatrixD gamma (fNa, 1);
        TMatrixD dgamma(fNa, 1);

        Double_t fWinitial;
        fWinitial = 0.0;
        SetInitialWeight(fWinitial);
        // for my example this fWinitial was not good; therefore :
        fWinitial = 1.0;

        Int_t ix;
        Double_t xiter;

        //--------   start  scanning weights   --------------------------
        // if full == kFALSE   only quantities necessary for the Gauss-Newton
        //                     iteration are calculated in SchmellCore
        // if full == kTRUE    in addition the unfolded distribution,
        //                     its covariance matrix and quantities like
        //                     Chisq, SpurAR, etc. are computed in SchmellCore
        //Bool_t full;
        //full = kFALSE;
        Int_t full;

        cout << "Schmelling : start loop over weights" << endl;

        dga2 = 1.e20;
        for (ix=0; ix<Nix; ix++)
        {
            xiter = pow(10.0,log10(xmin)+ix*dlogx) * fWinitial;

            //----------   start Gauss-Newton iteration   ----------------------
            numGiteration = 0;

            // if there was no convergence and the starting gamma was != 0
            // redo unfolding for the same weight starting with gamma = 0
            //
            Int_t gamma0 = 0;
            while (1)
	    {
              if (dga2 > EpsG)
	      {
                gamma0 = 1;
                gamma.Zero();
	      }

              dga2 = 1.e20;

              while (1)
              {
                dga2old = dga2;

                full = 0;
                SchmellCore(full, xiter, gamma, dgamma, dga2);

                gamma += dgamma;

                //cout << "Schmelling : ix, numGiteration, dga2, dga2old = "
                //     << ix << ",  " << numGiteration << ",  "
                //     << dga2 << ",  " << dga2old << endl;

                numGiteration += 1;

                // convergence
                if (dga2 < EpsG)
                    break;

                // no convergence
                if (numGiteration > MaxGiteration)
                    break;

                // gamma doesn't seem to change any more
                if (fabs(dga2-dga2old) < EpsG/100.)
                    break;
              }
              //----------   end Gauss-Newton iteration   ------------------------
              if (dga2<EpsG || gamma0 != 0) break;
	    }

            // if Gauss-Newton iteration has not converged
            // go to next weight
            if (dga2 > EpsG)
            {
                cout << "Schmelling : Gauss-Newton iteration has not converged;"
                    << "   numGiteration = " << numGiteration << endl;
                cout << "             ix, dga2, dga2old = " << ix << ",  "
                    << dga2 << ",  " << dga2old << endl;
                continue;
            }

            //cout << "Schmelling : Gauss-Newton iteration has converged;" << endl;
            //cout << "==================================================" << endl;
            //cout << "             numGiteration = " << numGiteration << endl;
            //cout << "             ix, dga2 = " << ix << ",  " << dga2 << endl;

            // calculated quantities which will be useful for determining
            // the best weight (Chisq, SpurAR, ...)
            //full = kTRUE;
            full = 1;
            SchmellCore(full, xiter, gamma, dgamma, dga2);

            // calculate difference between ideal and unfolded distribution
            Double_t D2bar = 0.0;
            for (UInt_t i = 0; i<fNb; i++)
            {
                Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
                D2bar += temp*temp;
            }

            SpAR(ix)  = SpurAR;
            SpSig(ix) = SpurSigma;
            chisq(ix) = Chisq;
            SecDer(ix) = SecDeriv;
            ZerDer(ix) = ZerDeriv;
            Entrop(ix) = Entropy;
            DAR2(ix)   = DiffAR2;
            Dsqbar(ix) = D2bar;

        }
        //----------   end of scanning weights   -------------------------------
        cout << "Schmelling : end of loop over weights" << endl;
        // plots ------------------------------
        for (ix=0; ix<Nix; ix++)
        {
            Double_t xbin = log10(xmin)+ix*dlogx;
            xiter = pow(10.0,xbin) * fWinitial;

            Int_t bin;
            bin = hBchisq->FindBin( xbin );
            hBchisq->SetBinContent(bin,chisq(ix));
            hBSpAR->SetBinContent(bin,SpAR(ix));
            hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
            hBSecDeriv->SetBinContent(bin,SecDer(ix));
            hBZerDeriv->SetBinContent(bin,ZerDer(ix));
            hBEntropy->SetBinContent(bin,Entrop(ix));
            hBDAR2->SetBinContent(bin,DAR2(ix));
            hBD2bar->SetBinContent(bin,Dsqbar(ix));

            if (ix > 0)
            {
                Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
                hBDSpAR->SetBinContent(bin,DSpAR);

                Double_t diff = SpSig(ix) - SpSig(ix-1);
                Double_t DSpSig = diff;
                hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);

                Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
                hBDEntropy->SetBinContent(bin,DEntrop);

                Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
                hBDSecDeriv->SetBinContent(bin,DSecDer);

                Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
                hBDZerDeriv->SetBinContent(bin,DZerDer);
            }
        }

        // Select best weight
        SelectBestWeight();

        if (ixbest < 0.0)
        {
            cout << "Schmelling : no solution found; " << endl;
            return kFALSE;
        }

        // do the unfolding using the best weight
        //full = kTRUE;


        xiter = pow(10.0,log10(xmin)+ixbest*dlogx) * fWinitial;

        //----------   start Gauss-Newton iteration   ----------------------
        numGiteration = 0;
        gamma.Zero();
        dga2 = 1.e20;

        while (1)
        {
            full = 1;
            SchmellCore(full, xiter, gamma, dgamma, dga2);
            gamma += dgamma;

            //cout << "Schmelling : sum(dgamma^2) = " << dga2 << endl;

            numGiteration += 1;

            if (numGiteration > MaxGiteration)
                break;

            if (dga2 < EpsG)
                break;
        }
        //----------   end Gauss-Newton iteration   ------------------------


        //-----------------------------------------------------------------------
        // termination stage
        // =================

        cout << "Schmelling : best solution found; " << endl;
        cout << "==================================" << endl;
        cout << "             xiter, ixbest, numGiteration, Chisq = "
            << xiter << ",  " << ixbest << ",  "
            << numGiteration << ",  " << Chisq << endl;

        //------------------------------------
        //..............................................
        // put unfolded distribution into fResult
        //     fResult(i,0)   value in bin i
        //     fResult(i,1)   error of value in bin i

        fNdf = SpurAR;
        fChisq = Chisq;

        for (UInt_t i=0; i<fNa; i++)
	{
          fChi2(i,0) = Chi2(i,0);
	}

        UInt_t iNdf   = (UInt_t) (fNdf+0.5);
        fProb = iNdf>0 ? TMath::Prob(fChisq, iNdf) : 0;

        fResult.ResizeTo(fNb, 5);
        for (UInt_t i=0; i<fNb; i++)
        {
            fResult(i, 0) = fVb(i,0);
            fResult(i, 1) = sqrt(fVbcov(i,i));
            fResult(i, 2) = 0.0;
            fResult(i, 3) = 0.0;
            fResult(i, 4) = 1.0; 
	}

        //--------------------------------------------------------

        cout << "Schmelling : gamma = " << endl;
        for (UInt_t i=0; i<fNa; i++)
            cout << gamma(i,0) << " \t";
        cout << endl;

        return kTRUE;
    }




    // -----------------------------------------------------------------------
    //
    // SchmellCore     main part of Schmellings calculations
    //
    void SchmellCore(Int_t full, Double_t &xiter, TMatrixD &gamma,
                     TMatrixD &dgamma, Double_t &dga2)
    {
        Double_t norm;
        TMatrixD wZdp_inv(fNa, fNa);
        TMatrixD d(fNb, 1);
        TMatrixD p(fNb, 1);

        //--------------------------------------------------------
        //--      determine the probability vector p


        TMatrixD v3(gamma, TMatrixD::kTransposeMult, fMigrat);
        TMatrixD v4(TMatrixD::kTransposed, v3);
        d = v4;
        Double_t dmax  = -1.e10;
        for (UInt_t j=0; j<fNb; j++)
            if (d(j,0)>dmax)
                dmax = d(j,0);

        Double_t psum = 0.0;
        for (UInt_t j=0; j<fNb; j++)
        {
            d(j,0) -= dmax;
            p(j,0)  = fVEps0(j)*exp(xiter*d(j,0));
            psum += p(j,0);
        }

        p *= 1.0/psum;

        //--      get the vector alpha

        TMatrixD alpha(fMigrat, TMatrixD::kMult, p);

        //--      determine the current normalization

        TMatrixD v2   (alpha, TMatrixD::kTransposeMult, fVacovInv);
        TMatrixD normb(v2,    TMatrixD::kMult, alpha);

        TMatrixD normc(v2,    TMatrixD::kMult, fVa);

        norm  = normc(0,0)/normb(0,0);

        //--------------------------------------------------------
        //--      determine the scaled slope vector s and Hessian H

        TMatrixD Zp(fNa,1);
        for (UInt_t i=0; i<fNa; i++)
        {
            Zp(i,0) = norm*alpha(i,0) - fVa(i,0);
            for (UInt_t k=0; k<fNa; k++)
                Zp(i,0) += gamma(k,0)*fVacov(k,i);
        }


        TMatrixD Q  (fNa, fNa);
        TMatrixD wZdp(fNa, fNa);
        for (UInt_t i=0; i<fNa; i++)
        {
            for (UInt_t j=0; j<fNa; j++)
            {
                Q(i,j) = - alpha(i,0)*alpha(j,0);
                for (UInt_t k=0; k<fNb; k++)
                    Q(i,j) += fMigrat(i,k)*fMigrat(j,k)*p(k,0);
                wZdp(i,j) = xiter*norm*Q(i,j) + fVacov(i,j);
            }
        }

        //--      invert H and calculate the next Newton step

        Double_t determ = 1.0;
        wZdp_inv = wZdp;
        wZdp_inv.Invert(&determ);

        if(determ == 0.0)
        {
            cout << "SchmellCore: matrix inversion for H failed" << endl;
            return;
        }


        dga2 = 0.0;
        for (UInt_t i=0; i<fNa; i++)
        {
            dgamma(i,0) = 0.0;
            for (UInt_t j=0; j<fNa; j++)
                dgamma(i,0) -= wZdp_inv(i,j)*Zp(j,0);
            dga2 += dgamma(i,0)*dgamma(i,0);
        }

        if (full == 0)
            return;

        //--------------------------------------------------------
        //--      determine chi2 and dNdf (#measurements ignored)
        Double_t dNdf = 0;
        for (UInt_t i=0; i<fNa; i++)
        {
            Chi2(i,0) = 0;
            for (UInt_t j=0; j<fNa; j++)
            {
                Chi2(i,0) += fVacov(i,j) * gamma(i,0) * gamma(j,0);
                dNdf       += fVacov(i,j) * wZdp_inv(j,i);
            }
        }
        Chisq = GetMatrixSumCol(Chi2, 0);
        SpurAR = fNa - dNdf;

        //-----------------------------------------------------
        // calculate the norm |AR - AR+|**2

        TMatrixD AR(fNa, fNa);
        DiffAR2 = 0.0;
        for (UInt_t i=0; i<fNa; i++)
        {
            for (UInt_t j=0; j<fNa; j++)
            {
                AR(i,j) = 0.0;
                for (UInt_t k=0; k<fNa; k++)
                    AR(i,j) +=  fVacov(i,k) * wZdp_inv(k,j);
                DiffAR2 += AR(i,j) * AR(i,j);
            }
        }

        //--------------------------------------------------------
        //--      fill distribution b(*)
        fVb = p;
        fVb *= norm;

        //--      determine the covariance matrix of b (very expensive)

        TMatrixD T(fNb,fNa);
        for (UInt_t i=0; i<fNb; i++)
        {
            for (UInt_t j=0; j<fNa; j++)
            {
                T(i,j) = 0.0;
                for (UInt_t k=0; k<fNa; k++)
                    T(i,j) += xiter*wZdp_inv(k,j)*(fMigrat(k,i)-alpha(k,0))*p(i,0);
            }
        }

        SpurSigma = CalcSpurSigma(T, norm);

        //--------------------------------------------------------

        //-----------------------------------------------------
        // calculate the second derivative squared

        SecDeriv = 0;
        for (UInt_t j=1; j<(fNb-1); j++)
        {
            Double_t temp =
                + 2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
                - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));
            SecDeriv += temp*temp;
        }

        ZerDeriv = 0;
        for (UInt_t j=0; j<fNb; j++)
            ZerDeriv += fVb(j,0) * fVb(j,0);

        //-----------------------------------------------------
        // calculate the entropy
        Entropy = 0;
        for (UInt_t j=0; j<fNb; j++)
            if (p(j,0) > 0.0)
                Entropy += p(j,0) * log( p(j,0) );

        //--------------------------------------------------------
    }


    // -----------------------------------------------------------------------
    //
    // Smooth migration matrix
    //              by fitting a function to the migration matrix
    //
    Bool_t SmoothMigrationMatrix(TH2D &hmigorig)
    {
        // copy histograms into matrices; the matrices will be used in fcnSmooth
        // ------------------------

      
      //cout << "MUnfold::SmoothMigrationMatrix : fNa, fNb = " << fNa << ",  " << fNb << endl;

      //cout << "MUnfold::SmoothMigrationMatrix:   fMigOrig = "  << endl;
      //cout << "========================================"  << endl;
        for (UInt_t i=0; i<fNa; i++)
        {
            for (UInt_t j=0; j<fNb; j++)
            {
                fMigOrig(i, j)     = hmigorig.GetBinContent(i+1, j+1);
	        //cout << fMigOrig(i, j) << " \t";
            }
            //cout << endl;
        }
      

        // ------------------------

      
        //cout << "MUnfold::SmoothMigrationMatrix :   fMigOrigerr2 = " << endl;
        //cout << "=============================================" << endl;
        for (UInt_t i=0; i<fNa; i++)
        {
            for (UInt_t j=0; j<fNb; j++)
            {
                fMigOrigerr2(i, j) =   hmigorig.GetBinError(i+1, j+1)
                    * hmigorig.GetBinError(i+1, j+1);

                //cout << fMigOrigerr2(i, j) << " \t";
            }
            //cout << endl;
        }
      

        // ------------------------
        // the number of free parameters (npar) is equal to 6:
        //     a0mean, a1mean, a2mean     
        //     <log10(Eest)>    = a0 + a1*log10(Etrue) + a2*SQR(log10(Etrue))
        //                                                     + log10(Etrue)  
        //     b0RMS,  b1RMS, b2RMS      
        //     RMS(log10(Eest)) = b0 + b1*log10(Etrue) + b2*SQR(log10(Etrue))
        // 
        UInt_t npar = 6;

        if (npar > 20)
        {
            cout << "MUnfold::SmoothMigrationMatrix : too many parameters,  npar = "
                << npar << endl;
            return kFALSE;
        }


        //..............................................
        // Find reasonable starting values for a0, a1 and b0, b1

        Double_t xbar   = 0.0;
        Double_t xxbar  = 0.0;

        Double_t ybarm  = 0.0;
        Double_t xybarm = 0.0;

        Double_t ybarR  = 0.0;
        Double_t xybarR = 0.0;

        Double_t Sum = 0.0;
        for (UInt_t j=0; j<fNb; j++)
        {
            Double_t x = (double)j + 0.5;

            Double_t meany = 0.0;
            Double_t RMSy  = 0.0;
            Double_t sum   = 0.0;
            for (UInt_t i=0; i<fNa; i++)
            {
                Double_t y = (double)i + 0.5;
                meany +=   y * fMigOrig(i, j);
                RMSy  += y*y * fMigOrig(i, j);
                sum   +=       fMigOrig(i, j);
            }
            if (sum > 0.0)
            {
                meany  = meany / sum;
                RMSy   =  RMSy / sum  - meany*meany;
                RMSy = sqrt(RMSy);

                Sum += sum;
                xbar  +=   x * sum;
                xxbar += x*x * sum;

                ybarm  +=   meany * sum;
                xybarm += x*meany * sum;

                ybarR  +=   RMSy  * sum;
                xybarR += x*RMSy  * sum;
            }
        }

        if (Sum > 0.0)
        {
            xbar   /= Sum;
            xxbar  /= Sum;

            ybarm  /= Sum;
            xybarm /= Sum;

            ybarR  /= Sum;
            xybarR /= Sum;
        }

        Double_t a1start = (xybarm - xbar*ybarm) / (xxbar - xbar*xbar);
        Double_t a0start = ybarm - a1start*xbar;
        a1start = a1start - 1.0;

        Double_t b1start = (xybarR - xbar*ybarR) / (xxbar - xbar*xbar);
        Double_t b0start = ybarR - b1start*xbar;

        cout << "MUnfold::SmoothMigrationMatrix : " << endl;
        cout << "============================" << endl;
        cout << "a0start, a1start = " << a0start << ",  " << a1start << endl;
        cout << "b0start, b1start = " << b0start << ",  " << b1start << endl;

        //..............................................
        // Set starting values and step sizes for parameters

        char name[20][100];
        Double_t vinit[20];
        Double_t  step[20];
        Double_t limlo[20];
        Double_t limup[20];
        Int_t    fix[20];

        sprintf(&name[0][0], "a0mean");
        vinit[0] = a0start;
        //vinit[0] = 1.0;
        step[0]  = 0.1;
        limlo[0] = 0.0;
        limup[0] = 0.0;
        fix[0]   = 0;

        sprintf(&name[1][0], "a1mean");
        vinit[1] = a1start;
        //vinit[1] = 0.0;
        step[1]  = 0.1;
        limlo[1] = 0.0;
        limup[1] = 0.0;
        fix[1]   = 0;

        sprintf(&name[2][0], "a2mean");
        vinit[2] = 0.0;
        step[2]  = 0.1;
        limlo[2] = 0.0;
        limup[2] = 0.0;
        fix[2]   = 0;

        sprintf(&name[3][0], "b0RMS");
        vinit[3] = b0start;
          //vinit[3] = 0.8;
        step[3]  = 0.1;
        limlo[3] = 1.e-20;
        limup[3] = 10.0;
        fix[3]   = 0;

        sprintf(&name[4][0], "b1RMS");
        vinit[4] = b1start;
        //vinit[4] = 0.0;
        step[4]  = 0.1;
        limlo[4] = 0.0;
        limup[4] = 0.0;
        fix[4]   = 0;

        sprintf(&name[5][0], "b2RMS");
        vinit[5] = 0.0;
        step[5]  = 0.1;
        limlo[5] = 0.0;
        limup[5] = 0.0;
        fix[5]   = 0;


        if ( CallMinuit(fcnSmooth, npar, name, vinit,
                        step, limlo, limup, fix) )
        {

            // ------------------------
            // fMigrat is the migration matrix to be used in the unfolding;
            // fMigrat, as set by the constructor, is overwritten
            //          by the smoothed migration matrix

            for (UInt_t i=0; i<fNa; i++)
                for (UInt_t j=0; j<fNb; j++)
                    fMigrat(i, j) = fMigSmoo(i, j);

            // ------------------------

            for (UInt_t i=0; i<fNa; i++)
                for (UInt_t j=0; j<fNb; j++)
                    fMigraterr2(i, j) = fMigSmooerr2(i, j);


            // normalize
            for (UInt_t j=0; j<fNb; j++)
            {
                Double_t sum = 0.0;
                for (UInt_t i=0; i<fNa; i++)
                    sum += fMigrat(i, j);

                //cout << "SmoothMigrationMatrix : normalization fMigrat; j, sum + "
                //     << j << ",  " << sum << endl;

                if (sum == 0.0)
                    continue;

                for (UInt_t i=0; i<fNa; i++)
                {
                    fMigrat(i, j)     /= sum;
                    fMigraterr2(i, j) /= (sum*sum);
                }
            }

            cout << "MUnfold::SmoothMigrationMatrix :   fMigrat = "  << endl;
            cout << "========================================"  << endl;
            for (UInt_t i=0; i<fNa; i++)
            {
                for (UInt_t j=0; j<fNb; j++)
                    cout << fMigrat(i, j) << " \t";
                cout << endl;
            }

	    /*
            cout << "MUnfold::SmoothMigrationMatrix :   fMigraterr2 = "  << endl;
            cout << "============================================"  << endl;
            for (UInt_t i=0; i<fNa; i++)
            {
                for (UInt_t j=0; j<fNb; j++)
                    cout << fMigraterr2(i, j) << " \t";
                cout << endl;
            }
	    */

            // ------------------------

            return kTRUE;
        }

        return kFALSE;
    }

    // -----------------------------------------------------------------------
    //
    // Prepare the call to MINUIT for the minimization of the function
    //         f = chi2*w/2 + reg, where reg is the regularization term
    //         reg is the sum the squared 2nd derivatives
    //                        of the unfolded distribution
    //
    // the corresponding fcn routine is 'fcnTikhonov2'
    //
    Bool_t Tikhonov2(TH1D &hb0)
    {
        // the number of free parameters (npar) is equal to
        // the number of bins (fNb) of the unfolded distribution minus 1,
        // because of the constraint that the total number of events
        // is fixed
        UInt_t npar = fNb-1;

        if (npar > 20)
        {
            cout << "MUnfold::Tikhonov2 : too many parameters,  npar = "
                << npar << ",  fNb = " << fNb << endl;
            return kFALSE;
        }

        // copy ideal distribution
        
        for (UInt_t i=1; i<=fNb; i++)
        {
            fhb0->SetBinContent(i, hb0.GetBinContent(i));
            fhb0->SetBinError  (i, hb0.GetBinError(i));
        }
        

        //---   start w loop -----------------------------------

        cout << "Tikhonov2 : start loop over weights" << endl;

        Int_t ix;
        Double_t xiter;

        for (ix=0; ix<Nix; ix++)
        {
            fW = pow(10.0,log10(xmin)+ix*dlogx);

            //..............................................
            // Set starting values and step sizes for parameters

            char name[20][100];
            Double_t vinit[20];
            Double_t  step[20];
            Double_t limlo[20];
            Double_t limup[20];
            Int_t    fix[20];

            for (UInt_t i=0; i<npar; i++)
            {
                sprintf(&name[i][0], "p%d", i+1);
                vinit[i] = fVEps0(i);
                step[i]  = fVEps0(i)/10;

                // lower and upper limits  (limlo=limup=0: no limits)
                //limlo[i] = 1.e-20;
                limlo[i] = -1.0;
                limup[i] = 1.0;
                fix[i]   = 0;
            }

            // calculate solution for the weight fW
            // flag non-convergence by chisq(ix) = 0.0
            chisq(ix) = 0.0;
            if ( CallMinuit(fcnTikhonov2, npar, name, vinit,
                            step, limlo, limup, fix) )
            {
                // calculate difference between ideal and unfolded distribution
                Double_t D2bar = 0.0;
                for (UInt_t i = 0; i<fNb; i++)
                {
                    Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
                    D2bar += temp*temp;
                }

                SpAR(ix)   = SpurAR;
                SpSig(ix)  = SpurSigma;
                chisq(ix)  = Chisq;
                SecDer(ix) = SecDeriv;
                ZerDer(ix) = ZerDeriv;
                Entrop(ix) = Entropy;
                DAR2(ix)   = DiffAR2;
                Dsqbar(ix) = D2bar;
            }
        }
        cout << "Tikhonov2 : end of loop over weights" << endl;

        // plots ------------------------------
        for (ix=0; ix<Nix; ix++)
        {
            // test whether minimization has converged
            if (chisq(ix) != 0.0)
            {
                xiter = pow(10.0,log10(xmin)+ix*dlogx);

                Int_t bin;
                bin = hBchisq->FindBin( log10(xiter) );
                hBchisq->SetBinContent(bin,chisq(ix));

                //hBSpAR->SetBinContent(bin,SpAR(ix));
                hBSpAR->SetBinContent(bin,0.0);

                hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
                hBSecDeriv->SetBinContent(bin,SecDer(ix));
                hBZerDeriv->SetBinContent(bin,ZerDer(ix));
                hBEntropy->SetBinContent(bin,Entrop(ix));

                //hBDAR2->SetBinContent(bin,DAR2(ix));
                hBDAR2->SetBinContent(bin,0.0);

                hBD2bar->SetBinContent(bin,Dsqbar(ix));

                if (ix > 0)
                {
                    //Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
                    //hBDSpAR->SetBinContent(bin,DSpAR);

                    Double_t diff = SpSig(ix) - SpSig(ix-1);
                    Double_t DSpSig = diff;
                    hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);

                    Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
                    hBDEntropy->SetBinContent(bin,DEntrop);

                    Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
                    hBDSecDeriv->SetBinContent(bin,DSecDer);

                    Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
                    hBDZerDeriv->SetBinContent(bin,DZerDer);
                }
            }
        }


        //---   end w loop -----------------------------------

        // Select best weight
        SelectBestWeight();

        cout << " Tikhonov2 : after SelectBestWeight" << endl;

        if (ixbest < 0.0)
        {
            cout << "Tikhonov2 : no result found; " << endl;
            return kFALSE;
        }

        cout << "Tikhonov2 : best result found; " << endl;
        cout << "===============================" << endl;
        cout << "           ixbest = " << ixbest << endl;


        // do a final unfolding using the best weight

        fW = pow(10.0,log10(xmin)+ixbest*dlogx);

        //..............................................
        // Set starting values and step sizes for parameters

        char name[20][100];
        Double_t vinit[20];
        Double_t  step[20];
        Double_t limlo[20];
        Double_t limup[20];
        Int_t    fix[20];

        for (UInt_t i=0; i<npar; i++)
        {
            sprintf(&name[i][0], "p%d", i+1);
            vinit[i] = fVEps0(i);
            step[i]  = fVEps0(i)/10;

            // lower and upper limits  (limlo=limup=0: no limits)
            //limlo[i] = 1.e-20;
            limlo[i] = -1.0;
            limup[i] = 1.0;
            fix[i]   = 0;
        }

        // calculate solution for the best weight
        CallMinuit(fcnTikhonov2, npar, name, vinit,
                   step, limlo, limup, fix);


        cout << "Tikhonov : Values for best weight " << endl;
        cout << "==================================" << endl;
        cout << "fW, ixbest, Chisq, SpurAR, SpurSigma, SecDeriv, ZerDeriv, Entrop, DiffAR2, D2bar = " << endl;
        cout << "      " << fW << ",  " << ixbest << ",  "
            << Chisq << ",  " << SpurAR << ",  "
            << SpurSigma << ",  " << SecDeriv << ",  " << ZerDeriv << ",  "
            << Entropy << ",  " << DiffAR2 << ",  "
            << Dsqbar(ixbest) << endl;

        return kTRUE;

    }


    // -----------------------------------------------------------------------
    //
    // Bertero :
    //
    // the unfolded distribution is calculated iteratively;
    // the number of iterations after which the iteration is stopped
    //            corresponds to the 'weight' in other methods
    // a small number of iterations corresponds to strong regularization
    // a high number to no regularization
    //
    // see : M.Bertero, INFN/TC-88/2 (1988)
    //       V.B.Anykeyev et al., NIM A303 (1991) 350
    //
    Bool_t Bertero(TH1D &hb0)
    {
        // copy ideal distribution
        
        for (UInt_t i=1; i<=fNb; i++)
        {
            fhb0->SetBinContent(i, hb0.GetBinContent(i));
            fhb0->SetBinError  (i, hb0.GetBinError(i));
        }
        

        TMatrixD bold(fNb, 1);
        bold.Zero();

        //----------------------------------------------------------

        Double_t db2 = 1.e20;


        TMatrixD aminusaest(fNa, 1);

        //-------   scan number of iterations   -----------------

        cout << "Bertero : start loop over number of iterations" << endl;

        Int_t ix;

        for (ix=0; ix<Nix; ix++)
        {
            Double_t xiter = pow(10.0,log10(xmin)+ix*dlogx);

            // calculate solution for the iteration number xiter
            BertCore(xiter);

            // calculate difference between ideal and unfolded distribution
            Double_t D2bar = 0.0;
            for (UInt_t i = 0; i<fNb; i++)
            {
                Double_t temp = fVb(i,0)-hb0.GetBinContent(i+1,0);
                D2bar += temp*temp;
            }

            SpAR(ix)   = SpurAR;
            SpSig(ix)  = SpurSigma;
            chisq(ix)  = Chisq;
            SecDer(ix) = SecDeriv;
            ZerDer(ix) = ZerDeriv;
            Entrop(ix) = Entropy;
            DAR2(ix)   = DiffAR2;
            Dsqbar(ix) = D2bar;

            db2 = 0.0;
            for (UInt_t i = 0; i<fNb; i++)
            {
                Double_t temp = fVb(i,0)-bold(i,0);
                db2 += temp*temp;
            }
            bold = fVb;

            //if (db2 < Epsdb2) break;

        }

        cout << "Bertero : end of loop over number of iterations" << endl;

        // plots ------------------------------
        for (ix=0; ix<Nix; ix++)
        {
            Double_t xiter = pow(10.0,log10(xmin)+ix*dlogx);

            Int_t bin;
            bin = hBchisq->FindBin( log10(xiter) );
            hBchisq->SetBinContent(bin,chisq(ix));
            hBSpAR->SetBinContent(bin,SpAR(ix));
            hBSpSig->SetBinContent(bin,SpSig(ix)/fSpurVacov);
            hBSecDeriv->SetBinContent(bin,SecDer(ix));
            hBZerDeriv->SetBinContent(bin,ZerDer(ix));
            hBEntropy->SetBinContent(bin,Entrop(ix));
            hBDAR2->SetBinContent(bin,DAR2(ix));
            hBD2bar->SetBinContent(bin,Dsqbar(ix));

            if (ix > 0)
            {
                Double_t DSpAR = SpAR(ix) - SpAR(ix-1);
                hBDSpAR->SetBinContent(bin,DSpAR);

                Double_t diff = SpSig(ix) - SpSig(ix-1);
                Double_t DSpSig = diff;
                hBDSpSig->SetBinContent(bin, DSpSig/fSpurVacov);

                Double_t DEntrop = Entrop(ix) - Entrop(ix-1);
                hBDEntropy->SetBinContent(bin,DEntrop);

                Double_t DSecDer = SecDer(ix) - SecDer(ix-1);
                hBDSecDeriv->SetBinContent(bin,DSecDer);

                Double_t DZerDer = ZerDer(ix) - ZerDer(ix-1);
                hBDZerDeriv->SetBinContent(bin,DZerDer);
            }
        }
        //-------   end of scan of number of iterations   -----------------

        // Select best weight
        SelectBestWeight();


        if (ixbest < 0.0)
        {
            cout << "Bertero : weight iteration has NOT converged; " << endl;
            return kFALSE;
        }

        cout << "Bertero : weight iteration has converged; " << endl;
        cout << "            ixbest = " << ixbest << endl;


        // do a final unfolding using the best weight

        // calculate solution for the iteration number xiter
        Double_t xiter = pow(10.0,log10(xmin)+ixbest*dlogx);
        BertCore(xiter);

        cout << "Bertero : Values for best weight " << endl;
        cout << "=================================" << endl;
        cout << "fW, ixbest, Chisq, SpurAR, SpurSigma, SecDeriv, ZerDeriv, Entrop, DiffAR2, D2bar = " << endl;
        cout << "      " << fW << ",  " << ixbest << ",  "
            << Chisq << ",  " << SpurAR << ",  "
            << SpurSigma << ",  " << SecDeriv << ",  " << ZerDeriv << ",  "
            << Entropy << ",  " << DiffAR2 << ",  "
            << Dsqbar(ixbest) << endl;

        //----------------

        fNdf   = SpurAR;
        fChisq = Chisq;

        for (UInt_t i=0; i<fNa; i++)
	{
          fChi2(i,0) = Chi2(i,0);
	}

        UInt_t iNdf   = (UInt_t) (fNdf+0.5);
        fProb = iNdf>0 ? TMath::Prob(fChisq, iNdf) : 0;


        fResult.ResizeTo(fNb, 5);
        for (UInt_t i=0; i<fNb; i++)
        {
            fResult(i, 0) = fVb(i,0);
            fResult(i, 1) = sqrt(fVbcov(i,i));
            fResult(i, 2) = 0.0;
            fResult(i, 3) = 0.0;
            fResult(i, 4) = 1.0;
        }

        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // main part of Bertero's calculations
    //
    Bool_t BertCore(Double_t &xiter)
    {
        // ignore eigen values which are smaller than EpsLambda
        TMatrixD G_inv(fNa, fNa);
        TMatrixD Gtil_inv(fNa, fNa);
        TMatrixD atil(fNb, fNa);
        TMatrixD aminusaest(fNa, 1);

        G_inv.Zero();
        Gtil_inv.Zero();
        SpurAR = 0.0;

        // -----   loop over eigen values   ------------------
        // calculate the approximate inverse of the matrix G
        //cout << "flaml = " << endl;

        UInt_t flagstart = 2;
        Double_t flaml=0;

        for (UInt_t l=0; l<fNa; l++)
        {
            if (EigenValue(l) < EpsLambda)
                continue;

            switch (flagstart)
            {
            case 1 :
                // This is the expression for f(lambda) if the initial C^(0)
                // is chosen to be zero
                flaml = 1.0 - pow(1.0-tau*EigenValue(l),xiter);
                break;

            case 2 :
                // This is the expression for f(lambda) if the initial C^(0)
                // is chosen to be equal to the measured distribution
                flaml =           1.0 - pow(1.0-tau*EigenValue(l),xiter)
                    + EigenValue(l) * pow(1.0-tau*EigenValue(l),xiter);
                break;
            }

            //  cout << flaml << ",  ";

            for (UInt_t m=0; m<fNa; m++)
            {
                for (UInt_t n=0; n<fNa; n++)
                {
                    G_inv(m,n)    += 1.0  /EigenValue(l) * Eigen(m,l)*Eigen(n,l);
                    Gtil_inv(m,n) += flaml/EigenValue(l) * Eigen(m,l)*Eigen(n,l);
                }
            }
            SpurAR += flaml;
        }
        //cout << endl;


        //cout << "Gtil_inv =" << endl;
        //for (Int_t m=0; m<fNa; m++)
        //{
        //  for (Int_t n=0; n<fNa; n++)
        //  {
        //    cout << Gtil_inv(m,n) << ",  ";
        //  }
        //  cout << endl;
        //}

        //-----------------------------------------------------
        // calculate the unfolded distribution b
        TMatrixD v2(fMigrat, TMatrixD::kTransposeMult, Gtil_inv);
        atil = v2;
        TMatrixD v4(atil, TMatrixD::kMult, fVa);
        fVb = v4;

        //-----------------------------------------------------
        // calculate AR and AR+
        TMatrixD AR(v2, TMatrixD::kMult, fMigrat);

        TMatrixD v3(fMigrat, TMatrixD::kTransposeMult, G_inv);
        TMatrixD ARplus(v3, TMatrixD::kMult, fMigrat);


        //-----------------------------------------------------
        // calculate the norm |AR - AR+|**2

        DiffAR2 = 0.0;
        for (UInt_t j=0; j<fNb; j++)
        {
            for (UInt_t k=0; k<fNb; k++)
            {
                Double_t tempo = AR(j,k) - ARplus(j,k);
                DiffAR2       += tempo*tempo;
            }
        }

        //-----------------------------------------------------
        // calculate the second derivative squared

        SecDeriv = 0;
        for (UInt_t j=1; j<(fNb-1); j++)
        {
            // temp = ( 2.0*fVb(j,0)-fVb(j-1,0)-fVb(j+1,0) ) / ( 2.0*fVb(j,0) );
            Double_t temp =    2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
                - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));
            SecDeriv += temp*temp;
        }

        ZerDeriv = 0;
        for (UInt_t j=0; j<fNb; j++)
            ZerDeriv += fVb(j,0) * fVb(j,0);

        //-----------------------------------------------------
        // calculate the entropy

        Double_t sumb = 0.0;
        for (UInt_t j=0; j<fNb; j++)
            sumb += fVb(j,0);

        TMatrixD p(fNb,1);
        p = fVb;
        if (sumb > 0.0)
            p *= 1.0/sumb;

        Entropy = 0;
        for (UInt_t j=0; j<fNb; j++)
            if (p(j,0) > 0.0)
                Entropy += p(j,0) * log( p(j,0) );

        //-----------------------------------------------------

        TMatrixD Gb(fMigrat, TMatrixD::kMult, fVb);
        aminusaest = fVa;
        aminusaest -= Gb;
    
        TMatrixD v1(1,fNa);
        for (UInt_t i=0; i<fNa; i++)
        {
            v1(0,i) = 0.0;
            for (UInt_t j=0; j<fNa; j++)
                v1(0,i) += aminusaest(j,0) * fVacovInv(j,i) ;
        }

        //-----------------------------------------------------
        // calculate error matrix fVbcov of unfolded distribution
        SpurSigma = CalcSpurSigma(atil);

        //-----------------------------------------------------
        // calculate the chi squared
        for (UInt_t i = 0; i<fNa; i++)
            Chi2(i,0) = v1(0,i) * aminusaest(i,0);

        Chisq = GetMatrixSumCol(Chi2,0);
        return kTRUE;
    }


    // -----------------------------------------------------------------------
    //
    // Calculate the matrix G = M * M(transposed)
    //           and its eigen values and eigen vectors
    //
    Bool_t CalculateG()
    {
        // Calculate matrix G = M*M(transposed)     (M = migration matrix)
        //           the matrix Eigen of the eigen vectors of G
        //           the vector EigenValues of the eigen values of G
        //           the parameter tau = 1/lambda_max
        //
        TMatrixD v5(TMatrixD::kTransposed, fMigrat);
        //TMatrixD G(fMigrat, TMatrixD::kMult, v5);
        G.Mult(fMigrat, v5);

        Eigen = G.EigenVectors(EigenValue);

        RankG = 0.0;
        for (UInt_t l=0; l<fNa; l++)
        {
            if (EigenValue(l) < EpsLambda) continue;
            RankG += 1.0;
        }

        tau = 1.0 / EigenValue(0);

        //  cout << "eigen values : " << endl;
        //  for (Int_t i=0; i<fNa; i++)
        //  {
        //    cout << EigenValue(i) << ",  ";
        //  }
        //  cout << endl;

        //cout << "eigen vectors : " << endl;
        //for (Int_t i=0; i<fNa; i++)
        //{
        //  cout << "               vector " << i << endl;
        //  for (Int_t j=0; j<fNa; j++)
        //  {
        //    cout << Eigen(j,i) << ",  ";
        //  }
        //  cout << endl;
        //}
        //cout << endl;

        //cout << "G =" << endl;
        //for (Int_t m=0; m<fNa; m++)
        //{
        //  for (Int_t n=0; n<fNa; n++)
        //  {
        //    cout << G(m,n) << ",  ";
        //  }
        //  cout << endl;
        //}

        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // Select the best weight
    //
    Bool_t SelectBestWeight()
    {
        //-------------------------------
        // select 'best' weight according to some criterion

        Int_t ix;

        Double_t DiffSpSigmax = -1.e10;
        Int_t    ixDiffSpSigmax = -1;

        Double_t DiffSigpointsmin = 1.e10;
        Int_t    ixDiffSigpointsmin = -1;

        Double_t DiffRankGmin = 1.e10;
        Int_t    ixDiffRankGmin = -1;

        Double_t D2barmin = 1.e10;
        Int_t    ixD2barmin = -1;

        Double_t DiffSpSig1min = 1.e10;
        Int_t    ixDiffSpSig1min = -1;


        Int_t ixmax = -1;

        // first loop over all weights :
        //       find smallest chi2
        Double_t chisqmin = 1.e20;
        for (ix=0; ix<Nix; ix++)
        {
            // consider only weights for which
            //  - unfolding was successful
            if (chisq(ix) != 0.0)
            {
                if (chisq(ix) < chisqmin)
                    chisqmin = chisq(ix);
            }
        }
        Double_t chisq0 = chisqmin > fVapoints ? chisqmin : fVapoints/2.0;

        // second loop over all weights :
        //        consider only weights for which chisq(ix) < chisq0
        ixbest = -1;
        for (ix=0; ix<Nix; ix++)
        {
            if (chisq(ix) != 0.0 && chisq(ix) < 2.0*chisq0)
            {
                // ixmax = highest weight with successful unfolding
                //         (Least squares solution)
                ixmax = ix;

                SpurSigma = SpSig(ix);
                SpurAR    = SpAR(ix);
                Chisq    = chisq(ix);
                D2bar     = Dsqbar(ix);

                //----------------------------------
                // search weight where SpurSigma changes most
                //                               (as a function of the weight)
                if (ix > 0  &&  chisq(ix-1) != 0.0)
                {
                    Double_t diff = SpSig(ix) - SpSig(ix-1);
                    if (diff > DiffSpSigmax)
                    {
                        DiffSpSigmax   = diff;
                        ixDiffSpSigmax = ix;
                    }
                }

                //----------------------------------
                // search weight where Chisq is close
                //        to the number of significant measurements
                Double_t DiffSigpoints = fabs(Chisq-fVapoints);

                if (DiffSigpoints < DiffSigpointsmin)
                {
                    DiffSigpointsmin   = DiffSigpoints;
                    ixDiffSigpointsmin = ix;
                }

                //----------------------------------
                // search weight where Chisq is close
                //        to the rank of matrix G
                Double_t DiffRankG = fabs(Chisq-RankG);

                if (DiffRankG < DiffRankGmin)
                {
                    DiffRankGmin   = DiffRankG;
                    ixDiffRankGmin = ix;
                }

                //----------------------------------
                // search weight where SpurSigma is close to 1.0
                Double_t DiffSpSig1 = fabs(SpurSigma/fSpurVacov-1.0);

                if (DiffSpSig1 < DiffSpSig1min)
                {
                    DiffSpSig1min   = DiffSpSig1;
                    ixDiffSpSig1min = ix;
                }

                //----------------------------------
                // search weight where D2bar is minimal

                if (D2bar < D2barmin)
                {
                    D2barmin   = D2bar;
                    ixD2barmin = ix;
                }

                //----------------------------------
            }
        }


        // choose solution where increase of SpurSigma is biggest
        //if ( DiffSpSigmax > 0.0)
        //  ixbest = ixDiffSpSigmax;
        //else
        //  ixbest = ixDiffSigpointsmin;

        // choose Least Squares Solution
	//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
        // ixbest = ixmax;
	//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

        // choose weight where chi2 is close to the number of significant
        // measurements
        // ixbest = ixDiffSigpointsmin;

        // choose weight where chi2 is close to the rank of matrix G
        // ixbest = ixDiffRankGmin;

        // choose weight where chi2 is close to the rank of matrix G
	//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$
           ixbest = ixDiffSpSig1min;
	//$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$

        cout << "SelectBestWeight : ixDiffSpSigmax, DiffSpSigmax = "
            << ixDiffSpSigmax << ",  " << DiffSpSigmax << endl;
        cout << "================== ixDiffSigpointsmin, DiffSigpointsmin = "
            << ixDiffSigpointsmin << ",  " << DiffSigpointsmin << endl;

        cout << "                   ixDiffRankGmin, DiffRankGmin = "
            << ixDiffRankGmin << ",  " << DiffRankGmin << endl;

        cout << "                   ixDiffSpSig1min, DiffSpSig1min = "
            << ixDiffSpSig1min << ",  " << DiffSpSig1min << endl;

        cout << "                   ixD2barmin, D2barmin = "
            << ixD2barmin << ",  " << D2barmin << endl;
        cout << "                   ixmax  = " << ixmax  << endl;
        cout << "                   ixbest = " << ixbest << endl;


        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // Draw the plots
    //
    Bool_t DrawPlots()
    {

        // in the plots, mark the weight which has been selected
        Double_t xbin = log10(xmin)+ixbest*dlogx;

        TMarker *m = new TMarker();
        m->SetMarkerSize(1);
        m->SetMarkerStyle(20);

        //-------------------------------------
        // draw the iteration plots
        TString ctitle = bintitle;
        ctitle += "Plots versus weight";
        TCanvas *c = new TCanvas("iter", ctitle, 900, 600);
        c->Divide(3,2);

        c->cd(1);
        hBchisq->Draw();
        gPad->SetLogy();
        hBchisq->SetXTitle("log10(iteration number)");
        hBchisq->SetYTitle("chisq");
        m->DrawMarker(xbin, log10(chisq(ixbest)));

        c->cd(2);
        hBD2bar->Draw();
        gPad->SetLogy();
        hBD2bar->SetXTitle("log10(iteration number)");
        hBD2bar->SetYTitle("(b_unfolded-b_ideal)**2");
        m->DrawMarker(xbin, log10(Dsqbar(ixbest)));

        /*
         c->cd(3);
         hBDAR2->Draw();
         gPad->SetLogy();
         strgx = "log10(iteration number)";
         strgy = "norm(AR-AR+)";
         hBDAR2->SetXTitle(strgx);
         hBDAR2->SetYTitle(strgy);
         m->DrawMarker(xbin, log10(DAR2(ixbest)));
         */

        c->cd(3);
        hBSecDeriv->Draw();
        hBSecDeriv->SetXTitle("log10(iteration number)");
        hBSecDeriv->SetYTitle("Second Derivative squared");
        m->DrawMarker(xbin, SecDer(ixbest));

        /*
         c->cd(8);
         hBDSecDeriv->Draw();
         strgx = "log10(iteration number)";
         strgy = "Delta(Second Derivative squared)";
         hBDSecDeriv->SetXTitle(strgx);
         hBDSecDeriv->SetYTitle(strgy);
         */

        /*
         c->cd(4);
         hBZerDeriv->Draw();
         strgx = "log10(iteration number)";
         strgy = "Zero Derivative squared";
         hBZerDeriv->SetXTitle(strgx);
         hBZerDeriv->SetYTitle(strgy);
         m->DrawMarker(xbin, ZerDer(ixbest));
         */

        /*
         c->cd(5);
         hBDZerDeriv->Draw();
         strgx = "log10(iteration number)";
         strgy = "Delta(Zero Derivative squared)";
         hBDZerDeriv->SetXTitle(strgx);
         hBDZerDeriv->SetYTitle(strgy);
         */

        c->cd(4);
        hBSpAR->Draw();
        hBSpAR->SetXTitle("log10(iteration number)");
        hBSpAR->SetYTitle("SpurAR");
        m->DrawMarker(xbin, SpAR(ixbest));


        /*
         c->cd(11);
         hBDSpAR->Draw();
         strgx = "log10(iteration number)";
         strgy = "Delta(SpurAR)";
         hBDSpAR->SetXTitle(strgx);
         hBDSpAR->SetYTitle(strgy);
         */

        c->cd(5);
        hBSpSig->Draw();
        hBSpSig->SetXTitle("log10(iteration number)");
        hBSpSig->SetYTitle("SpurSig/SpurC");
        m->DrawMarker(xbin, SpSig(ixbest)/fSpurVacov);

        /*
         c->cd(14);
         hBDSpSig->Draw();
         strgx = "log10(iteration number)";
         strgy = "Delta(SpurSig/SpurC)";
         hBDSpSig->SetXTitle(strgx);
         hBDSpSig->SetYTitle(strgy);
         */

        c->cd(6);
        hBEntropy->Draw();
        hBEntropy->SetXTitle("log10(iteration number)");
        hBEntropy->SetYTitle("Entropy");
        m->DrawMarker(xbin, Entrop(ixbest));

        /*
         c->cd(17);
         hBDEntropy->Draw();
         strgx = "log10(iteration number)";
         strgy = "Delta(Entropy)";
         hBDEntropy->SetXTitle(strgx);
         hBDEntropy->SetYTitle(strgy);
         */

        //-------------------------------------

        for (UInt_t i=0; i<fNa; i++)
        {
            fha->SetBinContent(i+1, fVa(i, 0)         );
            fha->SetBinError  (i+1, sqrt(fVacov(i, i)));

            for (UInt_t j=0; j<fNb; j++)
            {
                fhmig->SetBinContent(i+1, j+1, fMigOrig(i, j)           );
                fhmig->SetBinError  (i+1, j+1, sqrt(fMigOrigerr2(i, j)) );

                shmig->SetBinContent(i+1, j+1, fMigrat(i, j)           );
                shmig->SetBinError  (i+1, j+1, sqrt(fMigraterr2(i, j)) );
                shmigChi2->SetBinContent(i+1, j+1, fMigChi2(i, j)   );
            }
        }

        //PrintTH2Content(*shmig);
        //PrintTH2Content(*shmigChi2);

        //-------------------------------------
        CopyCol(*hprior, fVEps);
        CopyCol(*hb,     fVb);
        for (UInt_t i=0; i<fNb; i++)
            hb->SetBinError(i+1, sqrt(fVbcov(i, i)));

        PrintTH1Content(*hb);
        PrintTH1Error(*hb);

        //..............................................
        for (UInt_t i=0; i<fNa; i++)
            hEigen->SetBinContent(i+1, EigenValue(i));

        //..............................................
        // draw the plots
        TString cctitle = bintitle;
        cctitle +=  "Unfolding input";
        TCanvas *cc = new TCanvas("input", cctitle, 900, 600);
        cc->Divide(3, 2);

        // distribution to be unfolded
        cc->cd(1);
        fha->Draw();
        gPad->SetLogy();
        fha->SetXTitle("log10(E-est/GeV)");
        fha->SetYTitle("Counts");

        // superimpose unfolded distribution
        // hb->Draw("*HSAME");

        // prior distribution
        cc->cd(2);
        hprior->Draw();
        gPad->SetLogy();
        hprior->SetXTitle("log10(E-true/GeV)");
        hprior->SetYTitle("Counts");

        // migration matrix
        cc->cd(3);
        fhmig->Draw("box");
        fhmig->SetXTitle("log10(E-est/GeV)");
        fhmig->SetYTitle("log10(E-true/GeV)");

        // smoothed migration matrix
        cc->cd(4);
        shmig->Draw("box");
        shmig->SetXTitle("log10(E-est/GeV)");
        shmig->SetYTitle("log10(E-true/GeV)");

        // chi2 contributions for smoothing
        cc->cd(5);
        shmigChi2->Draw("box");
        shmigChi2->SetXTitle("log10(E-est/GeV)");
        shmigChi2->SetYTitle("log10(E-true/GeV)");

        // Eigenvalues of matrix M*M(transposed)
        cc->cd(6);
        hEigen->Draw();
        hEigen->SetXTitle("l");
        hEigen->SetYTitle("Eigen values Lambda_l of M*M(transposed)");


       //..............................................
        // draw the results
        TString crtitle = bintitle;
        crtitle +=  "Unfolding results";
        TCanvas *cr = new TCanvas("results", crtitle, 600, 600);
        cr->Divide(2, 2);

        // unfolded distribution
        cr->cd(1);
        hb->Draw();
        gPad->SetLogy();
        hb->SetXTitle("log10(E-true/GeV)");
        hb->SetYTitle("Counts");

	
        // covariance matrix of unfolded distribution
        cr->cd(2);
        TH1 *hbcov=DrawMatrixClone(fVbcov, "lego");
        hbcov->SetBins(fNb, hb->GetBinLowEdge(1), hb->GetBinLowEdge(fNb+1),
                       fNb, hb->GetBinLowEdge(1), hb->GetBinLowEdge(fNb+1));

        hbcov->SetName("hbcov");
        hbcov->SetTitle("Error matrix of distribution hb");
        hbcov->SetXTitle("log10(E-true/GeV)");
        hbcov->SetYTitle("log10(E-true/GeV)");
       
	
        // chi2 contributions
        cr->cd(3);
        TH1 *hchi2=DrawMatrixColClone(fChi2);
        hchi2->SetBins(fNa, fha->GetBinLowEdge(1), fha->GetBinLowEdge(fNa+1));

        hchi2->SetName("Chi2");
        hchi2->SetTitle("chi2 contributions");
        hchi2->SetXTitle("log10(E-est/GeV)");
        hchi2->SetYTitle("Chisquared");
	
	
        // ideal distribution
        
        cr->cd(4);
        fhb0->Draw();
        gPad->SetLogy();
        fhb0->SetXTitle("log10(E-true/GeV)");
        fhb0->SetYTitle("Counts");
        

        // superimpose unfolded distribution
        hb->Draw("*Hsame");
	

        return kTRUE;
    }


    // -----------------------------------------------------------------------
    //
    // Interface to MINUIT
    //
    //
    Bool_t CallMinuit(
                      void (*fcnx)(Int_t &, Double_t *, Double_t &, Double_t *, Int_t),
                      UInt_t npar, char name[20][100],
                      Double_t vinit[20], Double_t step[20],
                      Double_t limlo[20], Double_t limup[20], Int_t fix[20])
    {
        //
        // Be carefull: This is not thread safe
        //
        UInt_t maxpar = 100;

        if (npar > maxpar)
        {
            cout << "MUnfold::CallMinuit : too many parameters,  npar = " << fNb
                << ",   maxpar = " << maxpar << endl;
            return kFALSE;
        }

        //..............................................
        // Set the maximum number of parameters
        TMinuit minuit(maxpar);


        //..............................................
        // Set the print level
        // -1   no output except SHOW comands
        //  0   minimum output
        //  1   normal output (default)
        //  2   additional ouput giving intermediate results
        //  3   maximum output, showing progress of minimizations
        //
        Int_t printLevel = -1;
        minuit.SetPrintLevel(printLevel);

        //..............................................       
        // Printout for warnings
        //    SET WAR      print warnings
        //    SET NOW      suppress warnings
        Int_t errWarn;
        Double_t tmpwar = 0;
        minuit.mnexcm("SET NOW", &tmpwar, 0, errWarn);

        //..............................................
        // Set the address of the minimization function
        minuit.SetFCN(fcnx);

        //..............................................
        // Set starting values and step sizes for parameters
        for (UInt_t i=0; i<npar; i++)
        {
            if (minuit.DefineParameter(i, &name[i][0], vinit[i], step[i],
                                       limlo[i], limup[i]))
            {
                cout << "MUnfold::CallMinuit: Error in defining parameter "
                    << name << endl;
                return kFALSE;
            }
        }

        //..............................................
        //Int_t NumPars = minuit.GetNumPars();
        //cout << "MUnfold::CallMinuit :  number of free parameters = "
        //     << NumPars << endl;

        //..............................................
        // Minimization
        minuit.SetObjectFit(this);

        //..............................................
        // Error definition :
        //
        //    for chisquare function :
        //      up = 1.0   means calculate 1-standard deviation error
        //         = 4.0   means calculate 2-standard deviation error
        //
        //    for log(likelihood) function :
        //      up = 0.5   means calculate 1-standard deviation error
        //         = 2.0   means calculate 2-standard deviation error
        Double_t up = 1.0;
        minuit.SetErrorDef(up);



        // Int_t errMigrad;
        // Double_t tmp = 0;
        // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);


        //..............................................
        // fix a parameter
        for (UInt_t i=0; i<npar; i++)
        {
            if (fix[i] > 0)
            {
                Int_t parNo = i;
                minuit.FixParameter(parNo);
            }
        }

        //..............................................
        // Set maximum number of iterations (default = 500)
        Int_t maxiter = 100000;
        minuit.SetMaxIterations(maxiter);

        //..............................................
        // minimization by the method of Migrad
        // Int_t errMigrad;
        // Double_t tmp = 0;
        // minuit.mnexcm("MIGRAD", &tmp, 0, errMigrad);

        //..............................................       
        // same minimization as by Migrad
        // but switches to the SIMPLEX method if MIGRAD fails to converge
        Int_t errMinimize;
        Double_t tmp = 0;
        minuit.mnexcm("MINIMIZE", &tmp, 0, errMinimize);

        //..............................................       
        // check quality of minimization
        // istat = 0   covariance matrix not calculated
        //         1   diagonal approximation only (not accurate)
        //         2   full matrix, but forced positive-definite
        //         3   full accurate covariance matrix 
        //             (indication of normal convergence)
        Double_t fmin, fedm, errdef;
        Int_t    npari, nparx, istat;
        minuit.mnstat(fmin, fedm, errdef, npari, nparx, istat);

        if (errMinimize || istat < 3)
        {
            cout << "MUnfold::CallMinuit : Minimization failed" << endl;
            cout << "       fmin = " << fmin   << ",   fedm = "  << fedm
                << ",   errdef = "  << errdef << ",   istat = " << istat
                << endl;
            return kFALSE;
        }

        //..............................................
        // Minos error analysis
        // minuit.mnmnos();

        //..............................................       
        // Print current status of minimization
        // if nkode = 0    only function value
        //            1    parameter values, errors, limits
        //            2    values, errors, step sizes, internal values
        //            3    values, errors, step sizes, 1st derivatives
        //            4    values, paraboloc errors, MINOS errors
  
        //Int_t nkode = 4;
        //minuit.mnprin(nkode, fmin);

        //..............................................       
        // call fcn with IFLAG = 3 (final calculation : calculate p(chi2))
        // iflag = 1   initial calculations only
        //         2   calculate 1st derivatives and function
        //         3   calculate function only
        //         4   calculate function + final calculations
        const char *command = "CALL";
        Double_t iflag = 3;
        Int_t errfcn3;
        minuit.mnexcm(command, &iflag, 1, errfcn3); 

        return kTRUE;
    }

    // -----------------------------------------------------------------------
    //
    // Return the unfolded distribution
    //
    TMatrixD &GetVb() { return fVb; }

    // -----------------------------------------------------------------------
    //
    // Return the covariance matrix of the unfolded distribution
    //
    TMatrixD &GetVbcov() { return fVbcov; }

    // -----------------------------------------------------------------------
    //
    // Return the unfolded distribution + various errors
    //
    TMatrixD &GetResult() { return fResult; }

    // -----------------------------------------------------------------------
    //
    // Return the chisquared contributions
    //
    TMatrixD &GetChi2() { return fChi2; }

    // -----------------------------------------------------------------------
    //
    // Return the total chisquared
    //
    Double_t &GetChisq() { return fChisq; }

    // -----------------------------------------------------------------------
    //
    // Return the number of degrees of freedom
    //
    Double_t &GetNdf() { return fNdf; }

    // -----------------------------------------------------------------------
    //
    // Return the chisquared probability
    //
    Double_t &GetProb() { return fProb; }

    // -----------------------------------------------------------------------
    //
    // Return the smoothed migration matrix
    //
    TMatrixD &GetMigSmoo() { return fMigSmoo; }

    // -----------------------------------------------------------------------
    //
    // Return the error2 of the smoothed migration matrix
    //
    TMatrixD &GetMigSmooerr2() { return fMigSmooerr2; }

    // -----------------------------------------------------------------------
    //
    // Return the chi2 contributions for the smoothing
    //
    TMatrixD &GetMigChi2() { return fMigChi2; }
};
// end of definition of class MUnfold
///////////////////////////////////////////////////


// -----------------------------------------------------------------------
//
// fcnSmooth     (used by SmoothMigrationMatrix)
//
// is called by MINUIT
// for given values of the parameters it calculates the function 
//                                               to be minimized
//
void fcnSmooth(Int_t &npar, Double_t *gin, Double_t &f, 
               Double_t *par, Int_t iflag)
{
    MUnfold &gUnfold = *(MUnfold*)gMinuit->GetObjectFit();

    Double_t a0 = par[0];
    Double_t a1 = par[1];
    Double_t a2 = par[2];

    Double_t b0 = par[3];
    Double_t b1 = par[4];
    Double_t b2 = par[5];

    // loop over bins of log10(E-true)
    Double_t chi2 = 0.0;
    Int_t npoints = 0;
    Double_t func[20];

    for (UInt_t j=0; j<gUnfold.fNb; j++)
    {
        Double_t yj   = ((double)j) + 0.5;
        Double_t mean = a0 + a1*yj + a2*yj*yj + yj;
        Double_t RMS  = b0 + b1*yj + b2*yj*yj;

        if (RMS <= 0.0)
        {
            chi2 = 1.e20;
            break;
        }

        // loop over bins of log10(E-est)

        //.......................................
        Double_t function;
        Double_t sum=0.0;
        for (UInt_t i=0; i<gUnfold.fNa; i++)
        {
            Double_t xlow = (double)i;
            Double_t xup  = xlow + 1.0;
            Double_t xl  = (xlow- mean) / RMS;
            Double_t xu  = (xup - mean) / RMS;
            function = (TMath::Freq(xu) - TMath::Freq(xl));

            //cout << "i, xl, xu, function = " <<  i <<  ",  "  << xl << ",  "
            //     << xu  << ",  " << function << endl;

            if (function < 1.e-10)
                function = 0.0;

            func[i] = function;
            sum += function;
        }

        //      cout << "mean, RMS = "  << mean << ",  " << RMS
        //     << ",   j , sum of function = " << j << ",  " << sum << endl;

        //.......................................

        for (UInt_t i=0; i<gUnfold.fNa; i++)
        {
            if (sum != 0.0)
                func[i] /= sum;

            gUnfold.fMigSmoo(i,j) = func[i];
            gUnfold.fMigChi2(i,j) = 0.0;

            // if relative error is greater than 30 % ignore the point

            if (gUnfold.fMigOrig(i,j)     != 0 &&
                gUnfold.fMigOrigerr2(i,j) != 0 &&
                func[i] != 0  )
            {
                if (gUnfold.fMigOrigerr2(i,j)/
                    (gUnfold.fMigOrig(i,j)*gUnfold.fMigOrig(i,j)) <= 0.09)
                {
                    gUnfold.fMigChi2(i,j) =   ( gUnfold.fMigOrig(i,j) - func[i] )
                        * ( gUnfold.fMigOrig(i,j) - func[i] )
                        /   gUnfold.fMigOrigerr2(i,j);
                    chi2 += gUnfold.fMigChi2(i,j);
                    npoints += 1;
                }
            }
        }
        //.......................................

    }
    f = chi2;

    //cout << "fcnSmooth : f = " << f << endl;

    //--------------------------------------------------------------------
    // final calculations
    if (iflag == 3)
    {
        Int_t     NDF = npoints - npar;
        Double_t prob = TMath::Prob(chi2, NDF);

        cout << "fcnSmooth : npoints, chi2, NDF, prob = " << npoints << ",  ";
        cout << chi2 << ",  " << NDF << ",  " << prob << endl;
        cout << "=======================================" << endl;
    }
}

// -----------------------------------------------------------------------
//
// fcnTikhonov2     (used by Tikhonov2)
//
// is called by MINUIT
// for given values of the parameters it calculates the function F
// the free parameters are the first (fNb-1) elements
//                     of the normalized unfolded distribution
//
void fcnTikhonov2(Int_t &npar, Double_t *gin, Double_t &f,
                  Double_t *par, Int_t iflag)
{
    MUnfold &gUnfold = *(MUnfold*)gMinuit->GetObjectFit();

    // (npar+1) is the number of bins of the unfolded distribuition (fNb)
    //  npar    is the number of free parameters                    (fNb-1)

    UInt_t npar1 = npar + 1;

    UInt_t fNa = gUnfold.fNa;
    UInt_t fNb = gUnfold.fNb;
    if (npar1 != fNb)
    {
        cout << "fcnTikhonov2 : inconsistency in number of parameters; npar, fNb = ";
        cout << npar << ",  " << fNb << endl;
        //return;
    }
    npar1 = fNb;

    TMatrixD p(npar1, 1);
    TMatrixD &fVb = gUnfold.fVb;

    // p is the normalized unfolded distribution
    // sum(p(i,0)) from i=0 to npar is equal to 1
    Double_t sum = 0.0;
    for (Int_t i=0; i<npar; i++)
    {
        p(i,0) = par[i];
        sum += par[i];
    }
    p(npar,0) = 1.0 - sum;


    // all p(i,0) have to be greater than zero
    for (UInt_t i=0; i<npar1; i++)
        if (p(i,0) <= 0.0)
        {
            f = 1.e20;
            return;
        }

    //.......................
    // take least squares result for the normaliztion
    TMatrixD alpha(gUnfold.fMigrat, TMatrixD::kMult, p);

    //TMatrixD v4   (gUnfold.fVa, TMatrixD::kTransposeMult,
    //                                 gUnfold.fVacovInv);
    //TMatrixD norma(v4,  TMatrixD::kMult, gUnfold.fVa);

    TMatrixD v5   (alpha, TMatrixD::kTransposeMult, gUnfold.fVacovInv);
    TMatrixD normb(v5,    TMatrixD::kMult, alpha);

    TMatrixD normc(v5,    TMatrixD::kMult, gUnfold.fVa);

    Double_t norm  = normc(0,0)/normb(0,0);

    //.......................

    // b is the unnormalized unfolded distribution
    // sum(b(i,0)) from i=0 to npar is equal to norm
    //                       (the total number of events)
    for (UInt_t i=0; i<npar1; i++)
        fVb(i,0) = p(i,0) * norm;

    TMatrixD Gb(gUnfold.fMigrat, TMatrixD::kMult, fVb);
    TMatrixD v3(fNa, 1);
    v3 = gUnfold.fVa;
    v3 -= Gb;

    TMatrixD v1(1,fNa);
    for (UInt_t i=0; i<fNa; i++)
    {
        v1(0,i) = 0;
        for (UInt_t j=0; j<fNa; j++)
            v1(0,i) += v3(j,0) * gUnfold.fVacovInv(j,i) ;
    }

    for (UInt_t i = 0; i<fNa; i++)
        gUnfold.Chi2(i,0) = v1(0,i) * v3(i,0);

    gUnfold.Chisq = GetMatrixSumCol(gUnfold.Chi2,0);

    //-----------------------------------------------------
    // calculate 2nd derivative squared
    // regularization term (second derivative squared)
    gUnfold.SecDeriv = 0;
    for (UInt_t j=1; j<(fNb-1); j++)
     {
         const Double_t temp =
             + 2.0*(fVb(j+1,0)-fVb(j,0)) / (fVb(j+1,0)+fVb(j,0))
             - 2.0*(fVb(j,0)-fVb(j-1,0)) / (fVb(j,0)+fVb(j-1,0));

         gUnfold.SecDeriv += temp*temp;
     }

    gUnfold.ZerDeriv = 0;
    for (UInt_t j=0; j<fNb; j++)
        gUnfold.ZerDeriv += fVb(j,0) * fVb(j,0);

    f = gUnfold.Chisq/2 * gUnfold.fW + gUnfold.SecDeriv;

    //cout << "F=" << f      << " \tSecDeriv=" << gUnfold.SecDeriv
    //     << " \tchi2="
    //	  << gUnfold.Chisq << " \tfW=" << gUnfold.fW << endl;

    //--------------------------------------------------------------------
    // final calculations
    if (iflag == 3)
    {
        //..............................................
        // calculate external error matrix of the fitted parameters 'val'
        // extend it with the covariances for y=1-sum(val)
        Double_t emat[20][20];
        Int_t    ndim = 20;
        gMinuit->mnemat(&emat[0][0], ndim);

        Double_t covv = 0;
        for (UInt_t i=0; i<(gUnfold.fNb-1); i++)
        {
            Double_t cov = 0;
            for (UInt_t k=0; k<(gUnfold.fNb-1); k++)
                cov += emat[i][k];

            emat[i][gUnfold.fNb-1] = -cov;
            emat[gUnfold.fNb-1][i] = -cov;

            covv += cov;
        }
        emat[gUnfold.fNb-1][gUnfold.fNb-1] = covv;

        for (UInt_t i=0; i<gUnfold.fNb; i++)
            for (UInt_t k=0; k<gUnfold.fNb; k++)
                gUnfold.fVbcov(i,k) = emat[i][k] *norm*norm;

        //-----------------------------------------------------
        //..............................................
        // put unfolded distribution into fResult
        //     fResult(i,0)   value in bin i
        //     fResult(i,1)   error of value in bin i

        gUnfold.fResult.ResizeTo(gUnfold.fNb, 5);

        Double_t sum = 0;
        for (UInt_t i=0; i<(gUnfold.fNb-1); i++)
        {
            Double_t val;
            Double_t err;
            if (!gMinuit->GetParameter(i, val, err))
            {
                cout << "Error getting parameter #" << i << endl;
                return;
            }

            Double_t eplus;
            Double_t eminus;
            Double_t eparab;
            Double_t gcc;
            gMinuit->mnerrs(i, eplus, eminus, eparab, gcc);

            gUnfold.fVb(i, 0)     = val    * norm;

            gUnfold.fResult(i, 0) = val    * norm;
            gUnfold.fResult(i, 1) = eparab * norm;
            gUnfold.fResult(i, 2) = eplus  * norm;
            gUnfold.fResult(i, 3) = eminus * norm;
            gUnfold.fResult(i, 4) = gcc;
            sum += val;
        }
        gUnfold.fVb(gUnfold.fNb-1, 0)     = (1.0-sum) * norm;

        gUnfold.fResult(gUnfold.fNb-1, 0) = (1.0-sum) * norm;
        gUnfold.fResult(gUnfold.fNb-1, 1) =
            sqrt(gUnfold.fVbcov(gUnfold.fNb-1,gUnfold.fNb-1));
        gUnfold.fResult(gUnfold.fNb-1, 2) = 0;
        gUnfold.fResult(gUnfold.fNb-1, 3) = 0;
        gUnfold.fResult(gUnfold.fNb-1, 4) = 1;
        //..............................................

        //-----------------------------------------------------
        // calculate 0th derivative squared
        gUnfold.ZerDeriv = 0;
        for (UInt_t j=0; j<fNb; j++)
            gUnfold.ZerDeriv += fVb(j,0) * fVb(j,0);

        //-----------------------------------------------------
        // calculate the entropy

        gUnfold.Entropy = 0;
        for (UInt_t j=0; j<gUnfold.fNb; j++)
            if (p(j,0) > 0)
                gUnfold.Entropy += p(j,0) * log( p(j,0) );


        //-----------------------------------------------------
        // calculate SpurSigma
        gUnfold.SpurSigma = 0.0;
        for (UInt_t m=0; m<fNb; m++)
            gUnfold.SpurSigma += gUnfold.fVbcov(m,m);
        // cout << "SpurSigma =" << SpurSigma << endl;

        //-----------------------------------------------------
        gUnfold.SpurAR  = 0;
        gUnfold.DiffAR2 = 0;

        //-----------------------------------------------------
        gUnfold.fNdf   = gUnfold.fNa;
        gUnfold.fChisq   = gUnfold.Chisq;

        for (UInt_t i=0; i<fNa; i++)
	{
          gUnfold.fChi2(i,0) = gUnfold.Chi2(i,0);
	}


        UInt_t iNdf = (UInt_t) (gUnfold.fNdf+0.5);

        //cout << "fcnTikhonov2 : fW, chisq (from fcnF) = "
        //     << gUnfold.fW << ",  " << gUnfold.fChisq << endl;

        gUnfold.fProb = iNdf>0 ? TMath::Prob(gUnfold.fChisq, iNdf) : 0;
    }
}


// ======================================================
//
// SteerUnfold
//
void SteerUnfold(TString bintitle,
                 TH1D &ha,     TH2D &hacov, TH2D &hmig,
                 TH2D &hmigor, TH1D &hb0,   TH1D *hpr,
                 TH1D &hb)
{
    // ha      is the distribution to be unfolded
    // hacov   is the covariance matrix of the distribution ha
    // hmig    is the migration matrix;
    //         it is used in the unfolding unless it is overwritten
    //         by SmoothMigrationMatrix by the smoothed migration matrix
    // hmigor  is the migration matrix to be smoothed;
    //         the smoothed migration matrix will be used in the unfolding
    // hpr     the prior distribution
    //         it is only used if SetPriorInput(*hpr) is called   
    // hb      unfolded distribution

    //..............................................       
    // create an MUnfold object;
    // fill histograms into vectors and matrices

    MUnfold unfold(ha, hacov, hmig);
    unfold.bintitle = bintitle;

    //..............................................       
    // smooth the migration matrix;
    //        the smoothed migration matrix will be used in the unfolding
    //        hmig is the original (unsmoothed) migration matrix

    unfold.SmoothMigrationMatrix(hmigor);

    //..............................................       
    // define prior distribution (has always to be defined) 
    // the alternatives are  :

    // 1  SetPriorConstant() :   isotropic distribution
    // 2  SetPriorPower(gamma) : dN/dE = E^{-gamma}
    // 3  SetPriorInput(*hpr):   the distribution *hpr is used 
    // 4  SetPriorRebin(*ha) :   use rebinned histogram ha 

    UInt_t flagprior = 4;
    cout << "SteerUnfold : flagprior = " << flagprior << endl;
    cout << "==========================="<< endl;

    Bool_t errorprior=kTRUE;
    switch (flagprior)
    {
    case 1:
        unfold.SetPriorConstant();
        break;
    case 2:
        errorprior = unfold.SetPriorPower(1.5);
        break;
    case 3:
        if (!hpr)
        {
            cout << "Error: No hpr!" << endl;
            return;
        }
        errorprior = unfold.SetPriorInput(*hpr);
        break;
    case 4:
        errorprior = unfold.SetPriorRebin(ha);
        break;
    }
    if (!errorprior)
    {
        cout << "MUnfold::SetPrior... : failed.  flagprior = " ;
        cout << flagprior << endl;
        return;
    }

    //..............................................       
    // calculate the matrix G = M * M(transposed)
    //           M being the migration matrix

    unfold.CalculateG();

    //..............................................       
    // call steering routine for the actual unfolding;
    // the alternatives are :

    // 1  Schmelling : minimize the function Z by Gauss-Newton iteration;
    //                 the parameters to be fitted are gamma(i) = lambda(i)/w;

    // 2  Tikhonov2 :  regularization term is sum of (2nd deriv.)**2 ;
    //                 minimization by using MINUIT;
    //                 the parameters to be fitted are
    //                 the bin contents of the unfolded distribution

    // 3  Bertero:     minimization by iteration
    //                 

    UInt_t flagunfold = 1;
    cout << "SteerUnfold : flagunfold = "  << flagunfold << endl;
    cout << "============================" << endl;



    switch (flagunfold)
    {
    case 1: // Schmelling
        cout << "" << endl;
        cout << "Unfolding algorithm : Schmelling" << endl;
        cout << "================================" << endl;
        if (!unfold.Schmelling(hb0))
            cout << "MUnfold::Schmelling : failed." << endl;
        break;

    case 2: // Tikhonov2
        cout << "" << endl;
        cout << "Unfolding algorithm :   Tikhonov" << endl;
        cout << "================================" << endl;
        if (!unfold.Tikhonov2(hb0))
            cout << "MUnfold::Tikhonov2 : failed." << endl;
        break;

    case 3: // Bertero
        cout << "" << endl;
        cout << "Unfolding algorithm :    Bertero" << endl;
        cout << "================================" << endl;
        if (!unfold.Bertero(hb0))
            cout << "MUnfold::Bertero : failed." << endl;
        break;
    }    


    //..............................................       
    // Print fResult
    unfold.PrintResults();


    //..............................................       
    // Draw the plots
    unfold.DrawPlots();

    //..............................................       
    // get unfolded distribution
    TMatrixD &Vb    = unfold.GetVb();
    TMatrixD &Vbcov = unfold.GetVbcov();

    UInt_t fNb = unfold.fNb;

    for (UInt_t a=0; a<fNb; a++)
    {
      hb.SetBinContent(a+1, Vb(a,0));
      hb.SetBinError(a+1, sqrt(Vbcov(a, a)) );
    }

}

//__________________________________________________________________________
////////////////////////////////////////////////////////////////////////////
//                                                                        //
// doUnfolding    (to be called in the analysis)                          //
//                                                                        //
// arguments :                                                            //
//                                                                        //
//   INPUT                                                                //
//      TH2D &tobeunfolded : no.of excess events and its error            //
//                           vs. (E-est, Theta)                           //
//      TH3D &migration    : migration matrix (E-est, E_true, Theta)      //
//                                                                        //
//   OUITPUT                                                              //
//      TH2D &unfolded     : no.of excess events and its error            //
//                           vs. (E-true, Theta)                          //
//                                                                        //
// calls SteerUnfold to do the unfolding                                  //
//                                                                        //
// The "Theta" axis is only used to loop over the bins of theta           //
// and to do the unfolding for each bin of theta. Instead of theta        //
// any other variable (or a dummy variable) may be used.                  //
//                                                                        //
//                                                                        //
////////////////////////////////////////////////////////////////////////////

void doUnfolding(TH2D &tobeunfolded, TH3D &migration, TH2D &unfolded)
{
  TAxis &taxis  = *tobeunfolded.GetYaxis();
  Int_t numybins = taxis.GetNbins();

  for (Int_t m=1; m<=numybins; m++)
  {
    TString bintitle = "Bin ";
    bintitle += m;
    bintitle += ": ";

    // -----------------------------------------
    // ha : distribution to be unfolded

    TH1D &ha = *tobeunfolded.ProjectionX("", m, m, "e");
    TString title = bintitle;
    title += "E-est distr. to be unfolded";
    ha.SetNameTitle("ha", title);
    TAxis &aaxis = *ha.GetXaxis();
    Int_t na = aaxis.GetNbins();
    Double_t alow = aaxis.GetBinLowEdge(1);
    Double_t aup  = aaxis.GetBinLowEdge(na+1);

    PrintTH1Content(ha);
    PrintTH1Error(ha);

    // -----------------------------------------
    // covariance matrix of the distribution ha

    title = bintitle;
    title +=  "Error matrix of distribution ha";
    TH2D hacov("hacov", title, na, alow, aup, na, alow, aup);
    //MH::SetBinning(&hacov, &aaxis, &aaxis); 

    Double_t errmin = 3.0;
    for (Int_t i=1; i<=na; i++)
    {
      for (Int_t j=1; j<=na; j++)
      {
        hacov.SetBinContent(i, j, 0.0);
      }
      const Double_t content = ha.GetBinContent(i);
      const Double_t error2  = (ha.GetBinError(i))*(ha.GetBinError(i));
      if (content <= errmin  &&  error2 < errmin) 
        hacov.SetBinContent(i, i, errmin);
      else
        hacov.SetBinContent(i, i, error2);
    }

    //PrintTH2Content(hacov);
    

    // -----------------------------------------
    // migration matrix :
    //           x corresponds to measured quantity
    //           y corresponds to true quantity
    TH2D &hmig = *(TH2D*)migration.Project3D("yxe");
    title = bintitle;
    title += "Migration Matrix";
    hmig.SetNameTitle("Migrat", title);

    TAxis &aaxismig = *hmig.GetXaxis();
    Int_t namig = aaxismig.GetNbins();

    if (na != namig)
    {
      cout << "doUnfolding : binnings are incompatible; na, namig = "
           << na << ",  " << namig << endl;
      return;
    }

    TAxis &baxismig = *hmig.GetYaxis();
    Int_t nbmig = baxismig.GetNbins();
    Double_t blow = baxismig.GetBinLowEdge(1);
    Double_t bup  = baxismig.GetBinLowEdge(nbmig+1);

    PrintTH2Content(hmig);
    //PrintTH2Error(hmig);


    // -----------------------------------------
    // dummy ideal distribution

    Int_t nb = nbmig;

    title = bintitle;
    title += "Dummy Ideal distribution";
    TH1D hb0("dummyhb0", title, nb, blow, bup);;
    //MH::SetBinning(&hb0, &baxismig); 
    hb0.Sumw2();

    for (Int_t k=1; k<=nb; k++)
    {
        hb0.SetBinContent(k, 1.0/nb);
        hb0.SetBinError  (k, 0.1/nb);
    }

    //PrintTH1Content(hb0);

    // -----------------------------------------
    // here the prior distribution can be defined for the call
    // to SetPriorInput(*hpr)

    title = bintitle;
    title += "Dummy Prior distribution";
    TH1D hpr("hpr", title, nb, blow, bup);
    //MH::SetBinning(&hpr, &baxismig); 
    
    for (Int_t k=1; k<=nb; k++)
        hpr.SetBinContent(k, 1.0/nb);

    //PrintTH1Content(hpr);


    // -----------------------------------------
    // unfolded distribution


    title = bintitle;
    title += "Unfolded distribution";
    TH1D hb("hb", title, nb, blow, bup);
    //MH::SetBinning(&hb, &baxismig); 

    // -----------------------------------------
    SteerUnfold(bintitle, ha, hacov, hmig, hmig, hb0, &hpr, hb);

    for (Int_t k=1; k<=nb; k++)
    {
      Double_t content = hb.GetBinContent(k);
      Double_t error   = hb.GetBinError(k);

      unfolded.SetBinContent(k, m, content);
      unfolded.SetBinError(k, m, error);
    }    

    delete &ha;
    delete &hmig;
  }

}
//========================================================================//


////////////////////////////////////////////////////////////////////////////
//                                                                        //
// Main program   (for testing purposes)                                  //
//                                                                        //
//                defines the ideal distribution          (hb0)           //
//                defines the migration matrix            (hMigrat)       //
//                defines the distribution to be unfolded (hVa)           //
//                                                                        //
//                calls doUnfolding                                       //
//                      to do the unfolding                               //
//                                                                        //
////////////////////////////////////////////////////////////////////////////
void fluxunfold()
{
  // -----------------------------------------
  // migration matrix :
  //           x corresponds to measured quantity
  //           y corresponds to true quantity

    const Int_t  na   =   13;
    //const Int_t  na   =   18;
    const Axis_t alow = 0.25;
    const Axis_t aup  = 3.50;

    const Int_t  nb   =   11;
    //const Int_t  nb   =   22;
    const Axis_t blow = 0.50;
    const Axis_t bup  = 3.25;

    const Int_t  nc   =     1;
    const Axis_t clow =   0.0;
    const Axis_t cup  =   1.0;

    Int_t m = 1;

    TH3D migration("migration", "Migration Matrix", 
                   na, alow, aup, nb, blow, bup, nc, clow, cup);
    migration.Sumw2();

    // parametrize migration matrix as
    //     <log10(Eest)>      = a0 + a1*log10(Etrue) + a2*log10(Etrue)**2 
    //                             + log10(Etrue) 
    //     RMS( log10(Eest) ) = b0 + b1*log10(Etrue) + b2*log10(Etrue)**2
    Double_t a0 = 0.0;
    Double_t a1 = 0.0;
    Double_t a2 = 0.0;

    Double_t b0 = 0.26;
    Double_t b1 =-0.054;
    Double_t b2 = 0.0;

    TF1 f2("f2", "gaus(0)", alow, aup);
    f2.SetParName(0, "ampl");
    f2.SetParName(1, "mean");
    f2.SetParName(2, "sigma");

    // loop over log10(Etrue) bins
    TAxis &yaxis = *migration.GetYaxis();
    for (Int_t j=1; j<=nb; j++)
    {
        Double_t yvalue = yaxis.GetBinCenter(j);

        const Double_t mean  = a0 + a1*yvalue + a2*yvalue*yvalue + yvalue;
        const Double_t sigma = b0 + b1*yvalue + b2*yvalue*yvalue;
        const Double_t ampl  = 1./ ( sigma*TMath::Sqrt(2.0*TMath::Pi()));

        // gaus(0) is a substitute for [0]*exp( -0.5*( (x-[1])/[2] )**2 )
        f2.SetParameter(0, ampl);
        f2.SetParameter(1, mean);
        f2.SetParameter(2, sigma);

        // fill temporary 1-dim histogram with the function
        // fill the histogram using
        //    - either FillRandom
        //    - or using Freq
        TH1D htemp("temp", "temp", na, alow, aup);
        htemp.Sumw2();

        for (Int_t k=0; k<1000000; k++)
            htemp.Fill(f2.GetRandom());

        // copy it into the migration matrix
        Double_t sum = 0;
        for (Int_t i=1; i<=na; i++)
        {
            const Stat_t content = htemp.GetBinContent(i);
            migration.SetBinContent(i, j, m, content);
            sum += content;
        }

        // normalize migration matrix
        if (sum==0)
            continue;

        for (Int_t i=1; i<=na; i++)
        {
            const Stat_t content = migration.GetBinContent(i,j,m);
            migration.SetBinContent(i,j,m,       content/sum);
            migration.SetBinError  (i,j,m, sqrt(content)/sum);
        }
    }

    //PrintTH3Content(migration);
    //PrintTH3Error(migration);

    // -----------------------------------------
    // ideal distribution

    TH1D hb0("hb0", "Ideal distribution", nb, blow, bup);
    hb0.Sumw2();

    // fill histogram with random numbers according to 
    // an exponential function dN/dE = E^{-gamma}
    //    or with y = log10(E), E = 10^y :             
    //                          dN/dy = ln10 * 10^{y*(1-gamma)}     
    TF1 f1("f1", "pow(10.0, x*(1.0-[0]))", blow, bup);
    f1.SetParName(0,"gamma");
    f1.SetParameter(0, 2.7);

    // ntimes is the number of entries
    for (Int_t k=0; k<10000; k++)
        hb0.Fill(f1.GetRandom());

    // introduce energy threshold at 50 GeV

    const Double_t  lgEth = 1.70;
    const Double_t dlgEth = 0.09;

    for (Int_t j=1; j<=nb; j++)
    {
        const Double_t lgE = hb0.GetBinCenter(j);
        const Double_t c   = hb0.GetBinContent(j);
        const Double_t dc  = hb0.GetBinError(j);
        const Double_t f   = 1.0 / (1.0 + exp( -(lgE-lgEth)/dlgEth ));

        hb0.SetBinContent(j, f* c);
        hb0.SetBinError  (j, f*dc);
    }

    //PrintTH1Content(hb0);

    // -----------------------------------------
    // generate distribution to be unfolded (ha)
    // by smearing the ideal distribution  (hb0)
    //
    TH2D tobeunfolded("tobeunfolded", "Distribution to be unfolded", 
                      na, alow, aup, nc, clow, cup);
    tobeunfolded.Sumw2();

    for (Int_t i=1; i<=na; i++)
    {
        Double_t cont = 0;
        for (Int_t j=1; j<=nb; j++)
            cont += migration.GetBinContent(i, j, m) * hb0.GetBinContent(j);

        tobeunfolded.SetBinContent(i, m, cont);
        tobeunfolded.SetBinError(i, m, sqrt(cont));
    }

    //PrintTH2Content(tobeunfolded);
    //PrintTH2Error(tobeunfolded);

    // -----------------------------------------
    // unfolded distribution

    TH2D unfolded("unfolded", "Unfolded distribution", 
                  nb, blow, bup, nc, clow, cup);
    unfolded.Sumw2();

    // -----------------------------------------
    doUnfolding(tobeunfolded, migration, unfolded);

}
//========================================================================//