class TUnfoldSys: public TUnfold

 set all pointers to zero

 arguments:
    hist_A:  matrix that describes the migrations
    histmap: mapping of the histogram axes to the unfolding output
    regmode: global regularisation mode
 data members initialized to something different from zero:
    fDA2, fDAcol
 initialize TUnfold

 add a correlated error source
    sysError: alternative matrix or matrix of absolute/relative shifts
    name: name of the error source
    histmap: mapping of the histogram axes to the unfolding output
    mode: format of the error source

 performs background subtraction
 fY = fYData - fBgrIn
 fVyy = fVyyData + fBgrErrUncorr^2 + fBgrErrCorr * fBgrErrCorr#
 fVyyinv = fVyy^(-1)

 Define the input data for subsequent calls to DoUnfold(Double_t)
  input:   input distribution with errors
  scaleBias:  scale factor applied to the bias
  oneOverZeroError: for bins with zero error, this number defines 1/error.
 Return value: number of bins with bad error
                 +10000*number of unconstrained output bins
         Note: return values>=10000 are fatal errors,
               for the given input, the unfolding can not be done!
 Calls the SetInput metghod of the base class, then renames the input
 vectors fY and fVyy, then performs the background subtraction
 Data members modified:
   fYData,fY,fVyyData,fVyy,fVyyinvData,fVyyinv
 and those modified by TUnfold::SetInput()
 and those modified by DoBackgroundSubtraction()
LM: WARNING: Coverity detects here a false USE_AFTER_FREE for fY and fVyy
 the objects are deleted but then re-created immediatly afterwards in
  TUnfold::SetInput

 Store background source
   bgr:    background distribution with uncorrelated errors
   name:   name of this background source
   scale:  scale factor applied to the background
   scaleError: error on scale factor (correlated error)

 Data members modified:
   fBgrIn,fBgrErrUncorrIn,fBgrErrCorrIn
 and those modified by DoBackgroundSubtraction()
 save background source

 initialize pointers and TMaps

 delete all data members

 clear all data members which depend on the unfolding results

 calculations required for syst.error
 data members modified
    fEmatUncorrX, fEmatUncorrAx, fDeltaCorrX, fDeltaCorrAx

 get output error contribution from statistical fluctuations in A
   ematrix: output error matrix histogram
   binMap: see method GetEmatrix()
   clearEmat: set kTRUE to clear the histogram prior to adding the errors
 data members modified:
   fVYAx, fESparse, fEAtV, fErrorAStat

 propagate uncorrelated systematic errors to a covariance matrix
   m0,m1 : coefficients (matrices) for propagating the errors

 the error matrix is calculated by standard error propagation, where the
 derivative of the result vector X wrt the matrix A is given by

  dX_k / dA_ij  =  M0_kj * Z0_i  - M1_ki * Z1_j

 where:
   the matrices M0 and M1 are arguments to this function
   the vectors Z0, Z1 : GetDXDAZ()

 The matrix A is calculated from a matrix B as

    A_ij = B_ij / sum_k B_kj

 where k runs over additional indices of B, not present in A.
 (underflow and overflow bins, used for efficiency corrections)

 define:   Norm_j = sum_k B_kj   (data member fSumOverY)

 the derivative of A wrt this input matrix B is given by:

   dA_ij / dB_kj = (  delta_ik - A_ij ) * 1/Norm_j

 The covariance matrix Vxx is:

   Vxx_mn  = sum_ijlk [   (dX_m / dA_ij) * (dA_ij / dB_kj) * DB_kj
                        * (dX_n / dA_lj) * (dA_lj / dB_kj)  ]

 where DB_kj is the error on B_kj squared
 Simplify the sum over k:

   sum_k [ (dA_ij / dB_kj) * DB_kj * (dA_lj / dB_kj) ]
      =  sum_k [  ( delta_ik - A_ij ) * 1/Norm_j * DB_kj *
                * ( delta_lk - A_lj ) * 1/Norm_j ]
      =  sum_k [ ( delta_ik*delta_lk - delta_ik*A_lj - delta_lk*A_ij
                  + A_ij * A_lj ) * DB_kj / Norm_j^2 ]

 introduce normalized errors:  Rsq_kj = DB_kj / Norm_j^2
 after summing over k:
   delta_ik*delta_lk*Rsq_kj  ->    delta_il*Rsq_ij
   delta_ik*A_lj*Rsq_kj      ->    A_lj*Rsq_ij
   delta_lk*A_ij*Rsq_kj      ->    A_ij*Rsq_lj
   A_ij*A_lj*Rsq_kj          ->    A_ij*A_lj*sum_k(Rsq_kj)

 introduce sum of normalized errors squared:   SRsq_j = sum_k(Rsq_kj)

 Note: Rsq_ij is stored as  fDAinRelSq     (excludes extra indices of B)
   and SRsq_j is stored as  fDAinColRelSq  (sum includes all indices of B)

  Vxx_nm = sum_ijl [ (dX_m / dA_ij) * (dX_n / dA_lj)
     (delta_il*Rsq_ij - A_lj*Rsq_ij - A_ij*Rsq_lj + A_ij*A_lj *SRsq_j) ]

  Vxx_nm =    sum_j [ F_mj * F_nj * SRsq_j
            - sum_j [ G_mj * F_nj ]
            - sum_j [ F_mj * G_nj ]
            + sum_ij [  (dX_m / dA_ij) * (dX_n / dA_lj) * Rsq_ij ]

 where:
    F_mj = sum_i [ (dX_m / dA_ij) * A_ij ]
    G_mj = sum_i [ (dX_m / dA_ij) * Rsq_ij ]

 In order to avoid explicitly calculating the 3-dimensional tensor
 (dX_m/dA_ij) the sums are evaluated further, using
    dX_k / dA_ij  =  M0_kj * Z0_i  - M1_ki * Z1_j

   F_mj = M0_mj * (A# Z0)_j - (M1 A)_mj Z1_j
   G_mj = M0_mj * (Rsq# Z0)_j - (M1 Rsq)_mj Z1_j

 and

   sum_ij [ (dX_m/dA_ij) * (dX_n/dA_ij) * Rsq_ij ] =
      sum_j [ M0_mj * M0_nj *  [ sum_i (Z0_i)^2 * Rsq_ij ] ]
    + sum_i [ M1_mi * M1_ni *  [ sum_j (Z1_j)^2 * Rsq_ij ] ]
    - sum_i [ M1_mi * H_ni + M1_ni * H_mi]
 where:
   H_mi = Z0_i * sum_j [ M0_mj * Z1_j * Rsq_ij ]

 collect all contributions:
   Vxx_nm = r0 -r1 -r2 +r3 +r4 -r5 -r6
      r0 = sum_j [ F_mj * F_nj * SRsq_j ]
      r1 = sum_j [ G_mj * F_nj ]
      r2 = sum_j [ F_mj * G_nj ]
      r3 = sum_j [ M0_mj * M0_nj *  [ sum_i (Z0_i)^2 * Rsq_ij ] ]
      r4 = sum_i [ M1_mi * M1_ni *  [ sum_j (Z1_j)^2 * Rsq_ij ] ]
      r5 = sum_i [ M1_mi * H_ni ]
      r6 = sum_i [ M1_ni * H_mi ]

 calculate contributions containing matrices F and G
 r0,r1,r2

 set uncertainty on tau

 calculate systematic shift from a given source
    ematrix: output
    source: name of the error source
    binMap: see method GetEmatrix()

 get correlated shift induced by a background source
   delta: output shift vector histogram
   source: name of background source
   binMap: see method GetEmatrix()
   see PrepareSysError()

 calculate systematic shift from tau variation
    ematrix: output
    binMap: see method GetEmatrix()

 calculate systematic shift from a given source
    ematrix: output
    source: name of the error source
    binMap: see method GetEmatrix()
    clearEmat: set kTRUE to clear the histogram prior to adding the errors

 calculate systematic shift from a given background scale error
    ematrix: output
    source: name of the error source
    binMap: see method GetEmatrix()
    clearEmat: set kTRUE to clear the histogram prior to adding the errors

 calculate error matrix from error in regularisation parameter
    ematrix: output
    binMap: see method GetEmatrix()
    clearEmat: set kTRUE to clear the histogram prior to adding the errors

 calculate error matrix from error in input vector alone
    ematrix: output
    binMap: see method GetEmatrix()
    clearEmat: set kTRUE to clear the histogram prior to adding the errors

 propagate error matrix vyy to the result
    vyy: error matrix on input data fY
    ematrix: output
    binMap: see method GetEmatrix()
    clearEmat: set kTRUE to clear the histogram prior to adding the errors

 get total error including statistical error
    ematrix: output
    binMap: see method GetEmatrix()

 calculate total chi**2 including systematic errors

 sum over bins of *delta, as defined in binMap,fXToHist
   hist_delta: histogram to return summed vector
   delta: vector to sum and remap