From $ROOTSYS/tutorials/math/testUnfold3.C

// Author: Stefan Schmitt
// DESY, 14.10.2008

//  Version 17.0, changeto use the class TUnfoldDensity
//
//  History:
//  Version 16.1, parallel to changes in TUnfold
//     Version 16.0, parallel to changes in TUnfold
//     Version 15, simple example including background subtraction

#include <TMath.h>
#include <TCanvas.h>
#include <TRandom3.h>
#include <TFitter.h>
#include <TF1.h>
#include <TStyle.h>
#include <TVector.h>
#include <TGraph.h>

#include "TUnfoldDensity.h"

using namespace std;

/*
  This file is part of TUnfold.

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

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

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

///////////////////////////////////////////////////////////////////////
//
// Simple Test program for the class TUnfoldDensity
//
// 1-dimensional unfolding with background subtraction
//
//  the goal is to unfold the underlying "true" distribution of a variable Pt
//
//  the reconstructed Pt is measured in 24 bins from 4 to 28
//  the generator-level Pt is unfolded into 10 bins from 6 to 26
//    plus underflow bin from 0 to 6
//    plus overflow bin above 26
//  there are two background sources
//       bgr1 and bgr2
//  the signal has a finite trigger efficiency at a threshold of 8 GeV
//
//  one type of systematic error is studied, where the signal parameters are
//  changed
//
//  Finally, the unfolding is compared to a "bin-by-bin" correction method
//
///////////////////////////////////////////////////////////////////////

TRandom *rnd=0;

Double_t GenerateEvent(const Double_t *parm,
                       const Double_t *triggerParm,
                       Double_t *intLumi,
                       Bool_t *triggerFlag,
                       Double_t *ptGen,Int_t *iType)
{
   // generate an event
   // input:
   //      parameters for the event generator
   // return value:
   //      reconstructed Pt
   // output to pointers:
   //      integrated luminosity
   //      several variables only accessible on generator level
   //
   // the parm array defines the physical parameters
   //  parm[0]: background source 1 fraction
   //  parm[1]: background source 2 fraction
   //  parm[2]: lower limit of generated Pt distribution
   //  parm[3]: upper limit of generated Pt distribution
   //  parm[4]: exponent for Pt distribution signal
   //  parm[5]: exponent for Pt distribution background 1
   //  parm[6]: exponent for Pt distribution background 2
   //  parm[7]: resolution parameter a goes with sqrt(Pt)
   //  parm[8]: resolution parameter b goes with Pt
   //  triggerParm[0]: trigger threshold turn-on position
   //  triggerParm[1]: trigger threshold turn-on width
   //  triggerParm[2]: trigger efficiency for high pt
   //
   // intLumi is advanced bu 1 for each *generated* event
   // for data, several events may be generated, until one passes the trigger
   //
   // some generator-level quantities are also returned:
   //   triggerFlag: whether the event passed the trigger threshold
   //   ptGen: the generated pt
   //   iType: which type of process was simulated
   //
   // the "triggerFlag" also has another meaning:
   //   if(triggerFlag==0)   only triggered events are returned
   //
   // Usage to generate data events
   //      ptObs=GenerateEvent(parm,triggerParm,0,0,0)
   //
   // Usage to generate MC events
   //      ptGen=GenerateEvent(parm,triggerParm,&triggerFlag,&ptGen,&iType);
   //
   Double_t ptObs;
   Bool_t isTriggered=kFALSE;
   do {
      Int_t itype;
      Double_t ptgen;
      Double_t f=rnd->Rndm();
      // decide whether this is background or signal
      itype=0;
      if(f<parm[0]) itype=1;
      else if(f<parm[0]+parm[1]) itype=2;
      // generate Pt according to distribution  pow(ptgen,a)
      // get exponent
      Double_t a=parm[4+itype];
      Double_t a1=a+1.0;
      Double_t t=rnd->Rndm();
      if(a1 == 0.0) {
         Double_t x0=TMath::Log(parm[2]);
         ptgen=TMath::Exp(t*(TMath::Log(parm[3])-x0)+x0);
      } else {
         Double_t x0=pow(parm[2],a1);
         ptgen=pow(t*(pow(parm[3],a1)-x0)+x0,1./a1);
      }
      if(iType) *iType=itype;
      if(ptGen) *ptGen=ptgen;

      // smearing in Pt with large asymmetric tail
      Double_t sigma=
         TMath::Sqrt(parm[7]*parm[7]*ptgen+parm[8]*parm[8]*ptgen*ptgen);
      ptObs=rnd->BreitWigner(ptgen,sigma);

      // decide whether event was triggered
      Double_t triggerProb =
         triggerParm[2]/(1.+TMath::Exp((triggerParm[0]-ptObs)/triggerParm[1]));
      isTriggered= rnd->Rndm()<triggerProb;
      (*intLumi) ++;
   } while((!triggerFlag) && (!isTriggered));
   // return trigger decision
   if(triggerFlag) *triggerFlag=isTriggered;
   return ptObs;
}

//int main(int argc, char *argv[])
void testUnfold3()
{
  // switch on histogram errors
  TH1::SetDefaultSumw2();

  // show fit result
  gStyle->SetOptFit(1111);

  // random generator
  rnd=new TRandom3();

  // data and MC luminosities
  Double_t const lumiData= 30000;
  Double_t const lumiMC  =1000000;

  //========================
  // Step 1: define binning, book histograms

  // reconstructed pt (fine binning)
  Int_t const nDet=24;
  Double_t const xminDet=4.0;
  Double_t const xmaxDet=28.0;

  // generated pt (coarse binning)
  Int_t const nGen=10;
  Double_t const xminGen= 6.0;
  Double_t const xmaxGen=26.0;

  //==================================================================
  // book histograms
  //  (1) unfolding input: binning scheme is fine for detector,coarse for gen
  //   histUnfoldInput : reconstructed data, binning for unfolding
  //   histUnfoldMatrix : generated vs reconstructed distribution
  //   histUnfoldBgr1 : background source1, as predicted by MC
  //   histUnfoldBgr2 : background source2, as predicted by MC
  //      for systematic studies
  //   histUnfoldMatrixSys : generated vs reconstructed with different signal
  //
  //  (2) histograms required for bin-by-bin method
  //   histDetDATAbbb : reconstructed data for bin-by-bin
  //   histDetMCbbb : reconstructed MC for bin-by-bin
  //   histDetMCBgr1bbb : reconstructed MC bgr1 for bin-by-bin
  //   histDetMCBgr2bbb : reconstructed MC bgr2 for bin-by-bin
  //   histDetMCBgrPtbbb : reconstructed MC bgr from low/high pt for bin-by-bin
  //   histGenMCbbb : generated MC truth
  //     for systematic studies
  //   histDetMCSysbbb : reconstructed with changed trigger
  //   histGenMCSysbbb : generated MC truth
  //
  //  (3) monitoring and control
  //   histGenData : data truth for bias tests
  //   histDetMC : MC prediction

  // (1) create histograms required for unfolding
  TH1D *histUnfoldInput=
     new TH1D("unfolding input rec",";ptrec",nDet,xminDet,xmaxDet);
  TH2D *histUnfoldMatrix=
     new TH2D("unfolding matrix",";ptgen;ptrec",
              nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
  TH1D *histUnfoldBgr1=
     new TH1D("unfolding bgr1 rec",";ptrec",nDet,xminDet,xmaxDet);
  TH1D *histUnfoldBgr2=
     new TH1D("unfolding bgr2 rec",";ptrec",nDet,xminDet,xmaxDet);
  TH2D *histUnfoldMatrixSys=
     new TH2D("unfolding matrix sys",";ptgen;ptrec",
              nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);

  // (2) histograms required for bin-by-bin
  TH1D *histBbbInput=
     new TH1D("bbb input rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbSignalRec=
     new TH1D("bbb signal rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbBgr1=
     new TH1D("bbb bgr1 rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbBgr2=
     new TH1D("bbb bgr2 rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbBgrPt=
     new TH1D("bbb bgrPt rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbSignalGen=
     new TH1D("bbb signal gen",";ptgen",nGen,xminGen,xmaxGen);
  TH1D *histBbbSignalRecSys=
     new TH1D("bbb signal sys rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbBgrPtSys=
     new TH1D("bbb bgrPt sys rec",";ptrec",nGen,xminGen,xmaxGen);
  TH1D *histBbbSignalGenSys=
     new TH1D("bbb signal sys gen",";ptgen",nGen,xminGen,xmaxGen);

  // (3) control histograms
  TH1D *histDataTruth=
     new TH1D("DATA truth gen",";ptgen",nGen,xminGen,xmaxGen);
  TH1D *histDetMC=
     new TH1D("MC prediction rec",";ptrec",nDet,xminDet,xmaxDet);
  // ==============================================================
  // Step 2: generate data distributions
  //
  // data parameters: in real life these are unknown
  static Double_t parm_DATA[]={
     0.05, // fraction of background 1 (on generator level)
     0.05, // fraction of background 2 (on generator level)
     3.5, // lower Pt cut (generator level)
     100.,// upper Pt cut (generator level)
     -3.0,// signal exponent
     0.1, // background 1 exponent
     -0.5, // background 2 exponent
     0.2, // energy resolution a term
     0.01, // energy resolution b term
  };
  static Double_t triggerParm_DATA[]={8.0, // threshold is 8 GeV
                               4.0, // width is 4 GeV
                               0.95 // high Pt efficiency os 95%
  };

  Double_t intLumi=0.0;
  while(intLumi<lumiData) {
     Int_t iTypeGen;
     Bool_t isTriggered;
     Double_t ptGen;
     Double_t ptObs=GenerateEvent(parm_DATA,triggerParm_DATA,&intLumi,
                                  &isTriggered,&ptGen,&iTypeGen);
     if(isTriggered) {
        // (1) histogram for unfolding
        histUnfoldInput->Fill(ptObs);

        // (2) histogram for bin-by-bin
        histBbbInput->Fill(ptObs);
     }
     // (3) monitoring
     if(iTypeGen==0) histDataTruth->Fill(ptGen);
  }

  // ==============================================================
  // Step 3: generate default MC distributions
  //
  // MC parameters
  // default settings
  static Double_t parm_MC[]={
     0.05, // fraction of background 1 (on generator level)
     0.05, // fraction of background 2 (on generator level)
     3.5, // lower Pt cut (generator level)
     100.,// upper Pt cut (generator level)
     -4.0,// signal exponent !!! steeper than in data
          //                     to illustrate bin-by-bin bias
     0.1, // background 1 exponent
     -0.5, // background 2 exponent
     0.2, // energy resolution a term
     0.01, // energy resolution b term
  };
  static Double_t triggerParm_MC[]={8.0, // threshold is 8 GeV
                             4.0, // width is 4 GeV
                             0.95 // high Pt efficiency is 95%
  };

  // weighting factor to accomodate for the different luminosity in data and MC
  Double_t lumiWeight = lumiData/lumiMC;
  intLumi=0.0;
  while(intLumi<lumiMC) {
     Int_t iTypeGen;
     Bool_t isTriggered;
     Double_t ptGen;
     Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MC,&intLumi,&isTriggered,
                                  &ptGen,&iTypeGen);
     if(!isTriggered) ptObs=0.0;

     // (1) distribution required for unfolding

     if(iTypeGen==0) {
        histUnfoldMatrix->Fill(ptGen,ptObs,lumiWeight);
     } else if(iTypeGen==1) {
        histUnfoldBgr1->Fill(ptObs,lumiWeight);
     } else if(iTypeGen==2) {
        histUnfoldBgr2->Fill(ptObs,lumiWeight);
     }

     // (2) distributions for bbb
     if(iTypeGen==0) {
        if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
           histBbbSignalRec->Fill(ptObs,lumiWeight);
        } else {
           histBbbBgrPt->Fill(ptObs,lumiWeight);
        }
        histBbbSignalGen->Fill(ptGen,lumiWeight);
     } else if(iTypeGen==1) {
        histBbbBgr1->Fill(ptObs,lumiWeight);
     } else if(iTypeGen==2) {
        histBbbBgr2->Fill(ptObs,lumiWeight);
     }

     // (3) control distribution
     histDetMC ->Fill(ptObs,lumiWeight);
  }

  // ==============================================================
  // Step 4: generate MC distributions for systematic study
  // example: study dependence on initial signal shape
  //   -> BGR distributions do not change
  static Double_t parm_MC_SYS[]={
     0.05, // fraction of background: unchanged
     0.05, // fraction of background: unchanged
     3.5, // lower Pt cut (generator level)
     100.,// upper Pt cut (generator level)
     -2.0, // signal exponent changed
     0.1, // background 1 exponent
     -0.5, // background 2 exponent
     0.2, // energy resolution a term
     0.01, // energy resolution b term
  };

  intLumi=0.0;
  while(intLumi<lumiMC) {
     Int_t iTypeGen;
     Bool_t isTriggered;
     Double_t ptGen;
     Double_t ptObs=GenerateEvent(parm_MC_SYS,triggerParm_MC,&intLumi,
                                  &isTriggered,&ptGen,&iTypeGen);
     if(!isTriggered) ptObs=0.0;

     // (1) for unfolding
     if(iTypeGen==0) {
        histUnfoldMatrixSys->Fill(ptGen,ptObs);
     }

     // (2) for bin-by-bin
     if(iTypeGen==0) {
        if((ptGen>=xminGen)&&(ptGen<xmaxGen)) {
           histBbbSignalRecSys->Fill(ptObs);
        } else {
           histBbbBgrPtSys->Fill(ptObs);
        }
        histBbbSignalGenSys->Fill(ptGen);
     }
  }

  // this method is new in version 16 of TUnfold
  cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";

  //========================
  // Step 5: unfolding
  //

  // here one has to decide about the regularisation scheme

  // regularize curvature
  TUnfold::ERegMode regMode =
     TUnfold::kRegModeCurvature;

 // preserve the area
  TUnfold::EConstraint constraintMode=
     TUnfold::kEConstraintArea;

  // bin content is divided by the bin width
  TUnfoldDensity::EDensityMode densityFlags=
     TUnfoldDensity::kDensityModeBinWidth;

  // set up matrix of migrations
  TUnfoldDensity unfold(histUnfoldMatrix,TUnfold::kHistMapOutputHoriz,
                        regMode,constraintMode,densityFlags);

  // define the input vector (the measured data distribution)
  unfold.SetInput(histUnfoldInput);

  // subtract background, normalized to data luminosity
  //  with 10% scale error each
  Double_t scale_bgr=1.0;
  Double_t dscale_bgr=0.1;
  unfold.SubtractBackground(histUnfoldBgr1,"background1",scale_bgr,dscale_bgr);
  unfold.SubtractBackground(histUnfoldBgr2,"background2",scale_bgr,dscale_bgr);

  // add systematic error
  unfold.AddSysError(histUnfoldMatrixSys,"signalshape_SYS",
                     TUnfold::kHistMapOutputHoriz,
                     TUnfoldSys::kSysErrModeMatrix);

  // run the unfolding
  Int_t nScan=30;
  TSpline *logTauX,*logTauY;
  TGraph *lCurve;
  // this method scans the parameter tau and finds the kink in the L curve
  // finally, the unfolding is done for the best choice of tau
  Int_t iBest=unfold.ScanLcurve(nScan,0.,0.,&lCurve,&logTauX,&logTauY);

  cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
      <<" / "<<unfold.GetNdf()<<"\n";

  // save graphs with one point to visualize best choice of tau
  Double_t t[1],x[1],y[1];
  logTauX->GetKnot(iBest,t[0],x[0]);
  logTauY->GetKnot(iBest,t[0],y[0]);
  TGraph *bestLcurve=new TGraph(1,x,y);
  TGraph *bestLogTauLogChi2=new TGraph(1,t,x);


  //===========================
  // Step 6: retreive unfolding results

  // get unfolding output
  // includes the statistical and background errors
  // but not the other systematic uncertainties
  TH1 *histUnfoldOutput=unfold.GetOutput("PT(unfold,stat+bgrerr)");

  // retreive error matrix of statistical errors
  TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");

  // retreive full error matrix
  // This includes all systematic errors
  TH2 *histEmatTotal=unfold.GetEmatrixTotal("unfolding total error matrix");

  // create two copies of the unfolded data, one with statistical errors
  // the other with total errors
  TH1 *histUnfoldStat=new TH1D("PT(unfold,staterr)",";Pt(gen)",
                                nGen,xminGen,xmaxGen);
  TH1 *histUnfoldTotal=new TH1D("PT(unfold,totalerr)",";Pt(gen)",
                                 nGen,xminGen,xmaxGen);
  for(Int_t i=0;i<nGen+2;i++) {
     Double_t c=histUnfoldOutput->GetBinContent(i);

     // histogram with unfolded data and stat errors
     histUnfoldStat->SetBinContent(i,c);
     histUnfoldStat->SetBinError
        (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i)));

     // histogram with unfolded data and total errors
     histUnfoldTotal->SetBinContent(i,c);
     histUnfoldTotal->SetBinError
        (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i)));
  }

  // create histogram with correlation matrix
  TH2D *histCorr=new TH2D("Corr(total)",";Pt(gen);Pt(gen)",
                          nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
  for(Int_t i=0;i<nGen+2;i++) {
     Double_t ei,ej;
     ei=TMath::Sqrt(histEmatTotal->GetBinContent(i,i));
     if(ei<=0.0) continue;
     for(Int_t j=0;j<nGen+2;j++) {
        ej=TMath::Sqrt(histEmatTotal->GetBinContent(j,j));
        if(ej<=0.0) continue;
        histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej);
     }
  }

  // retreive bgr source 1
  TH1 *histDetNormBgr1=unfold.GetBackground("bgr1 normalized",
                                             "background1");
  TH1 *histDetNormBgrTotal=unfold.GetBackground("bgr total normalized");

  //========================
  // Step 7: plots

  TCanvas output;
  output.Divide(3,2);
  output.cd(1);
  // data, MC prediction, background
  histUnfoldInput->SetMinimum(0.0);
  histUnfoldInput->Draw("E");
  histDetMC->SetMinimum(0.0);
  histDetMC->SetLineColor(kBlue);
  histDetNormBgrTotal->SetLineColor(kRed);
  histDetNormBgr1->SetLineColor(kCyan);
  histDetMC->Draw("SAME HIST");
  histDetNormBgr1->Draw("SAME HIST");
  histDetNormBgrTotal->Draw("SAME HIST");

  output.cd(2);
  // unfolded data, data truth, MC truth
  histUnfoldTotal->SetMinimum(0.0);
  histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5);
  // outer error: total error
  histUnfoldTotal->Draw("E");
  // middle error: stat+bgr
  histUnfoldOutput->Draw("SAME E1");
  // inner error: stat only
  histUnfoldStat->Draw("SAME E1");
  histDataTruth->Draw("SAME HIST");
  histBbbSignalGen->SetLineColor(kBlue);
  histBbbSignalGen->Draw("SAME HIST");

  output.cd(3);
  // unfolding matrix
  histUnfoldMatrix->SetLineColor(kBlue);
  histUnfoldMatrix->Draw("BOX");

  // show tau as a function of chi**2
  output.cd(4);
  logTauX->Draw();
  bestLogTauLogChi2->SetMarkerColor(kRed);
  bestLogTauLogChi2->Draw("*");

  // show the L curve
  output.cd(5);
  lCurve->Draw("AL");
  bestLcurve->SetMarkerColor(kRed);
  bestLcurve->Draw("*");

  // show correlation matrix
  output.cd(6);
  histCorr->Draw("BOX");

  output.SaveAs("testUnfold3.ps");

  //============================================================
  // step 8: compare results to the so-called bin-by-bin "correction"

  std::cout<<"bin truth   result   (stat)   (bgr)    (sys)\n";
  std::cout<<"===+=====+=========+========+========+=======\n";
  for(Int_t i=1;i<=nGen;i++) {
     // data contribution in this bin
     Double_t data=histBbbInput->GetBinContent(i);
     Double_t errData=histBbbInput->GetBinError(i);

     // subtract background contributions
     Double_t data_bgr=data
        - scale_bgr*(histBbbBgr1->GetBinContent(i)
                     + histBbbBgr2->GetBinContent(i)
                     + histBbbBgrPt->GetBinContent(i));
     Double_t errData_bgr=TMath::Sqrt
        (errData*errData+
         TMath::Power(dscale_bgr*histBbbBgr1->GetBinContent(i),2)+
         TMath::Power(scale_bgr*histBbbBgr1->GetBinError(i),2)+
         TMath::Power(dscale_bgr*histBbbBgr2->GetBinContent(i),2)+
         TMath::Power(scale_bgr*histBbbBgr2->GetBinError(i),2)+
         TMath::Power(dscale_bgr*histBbbBgrPt->GetBinContent(i),2)+
         TMath::Power(scale_bgr*histBbbBgrPt->GetBinError(i),2));
     // "correct" the data, using the Monte Carlo and neglecting off-diagonals
     Double_t fCorr=(histBbbSignalGen->GetBinContent(i)/
                     histBbbSignalRec->GetBinContent(i));
     Double_t data_bbb= data_bgr *fCorr;
     // stat only error
     Double_t errData_stat_bbb = errData*fCorr;
     // stat plus background subtraction
     Double_t errData_statbgr_bbb = errData_bgr*fCorr;

     // estimate systematic error by repeating the exercise
     // using the MC with systematic shifts applied
     Double_t fCorr_sys=(histBbbSignalGenSys->GetBinContent(i)/
                         histBbbSignalRecSys->GetBinContent(i));
     Double_t shift_sys= data_bgr*fCorr_sys - data_bbb;

     // add systematic shift quadratically and get total error
     Double_t errData_total_bbb=
        TMath::Sqrt(errData_statbgr_bbb*errData_statbgr_bbb
                    +shift_sys*shift_sys);

     // get results from real unfolding
     Double_t data_unfold= histUnfoldStat->GetBinContent(i);
     Double_t errData_stat_unfold=histUnfoldStat->GetBinError(i);
     Double_t errData_statbgr_unfold=histUnfoldOutput->GetBinError(i);
     Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);

     // compare

     std::cout<<TString::Format
        ("%3d %5.0f %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (unfolding)",
         i,histDataTruth->GetBinContent(i),data_unfold,
         errData_stat_unfold,TMath::Sqrt(errData_statbgr_unfold*
                                         errData_statbgr_unfold-
                                         errData_stat_unfold*
                                         errData_stat_unfold),
         TMath::Sqrt(errData_total_unfold*
                     errData_total_unfold-
                     errData_statbgr_unfold*
                     errData_statbgr_unfold))<<"\n";
     std::cout<<TString::Format
        ("          %8.1f +/-%5.1f +/-%5.1f +/-%5.1f (bin-by-bin)",
         data_bbb,errData_stat_bbb,TMath::Sqrt(errData_statbgr_bbb*
                                               errData_statbgr_bbb-
                                               errData_stat_bbb*
                                               errData_stat_bbb),
         TMath::Sqrt(errData_total_bbb*
                     errData_total_bbb-
                     errData_statbgr_bbb*
                     errData_statbgr_bbb))
              <<"\n\n";
  }
}
 testUnfold3.C:1
 testUnfold3.C:2
 testUnfold3.C:3
 testUnfold3.C:4
 testUnfold3.C:5
 testUnfold3.C:6
 testUnfold3.C:7
 testUnfold3.C:8
 testUnfold3.C:9
 testUnfold3.C:10
 testUnfold3.C:11
 testUnfold3.C:12
 testUnfold3.C:13
 testUnfold3.C:14
 testUnfold3.C:15
 testUnfold3.C:16
 testUnfold3.C:17
 testUnfold3.C:18
 testUnfold3.C:19
 testUnfold3.C:20
 testUnfold3.C:21
 testUnfold3.C:22
 testUnfold3.C:23
 testUnfold3.C:24
 testUnfold3.C:25
 testUnfold3.C:26
 testUnfold3.C:27
 testUnfold3.C:28
 testUnfold3.C:29
 testUnfold3.C:30
 testUnfold3.C:31
 testUnfold3.C:32
 testUnfold3.C:33
 testUnfold3.C:34
 testUnfold3.C:35
 testUnfold3.C:36
 testUnfold3.C:37
 testUnfold3.C:38
 testUnfold3.C:39
 testUnfold3.C:40
 testUnfold3.C:41
 testUnfold3.C:42
 testUnfold3.C:43
 testUnfold3.C:44
 testUnfold3.C:45
 testUnfold3.C:46
 testUnfold3.C:47
 testUnfold3.C:48
 testUnfold3.C:49
 testUnfold3.C:50
 testUnfold3.C:51
 testUnfold3.C:52
 testUnfold3.C:53
 testUnfold3.C:54
 testUnfold3.C:55
 testUnfold3.C:56
 testUnfold3.C:57
 testUnfold3.C:58
 testUnfold3.C:59
 testUnfold3.C:60
 testUnfold3.C:61
 testUnfold3.C:62
 testUnfold3.C:63
 testUnfold3.C:64
 testUnfold3.C:65
 testUnfold3.C:66
 testUnfold3.C:67
 testUnfold3.C:68
 testUnfold3.C:69
 testUnfold3.C:70
 testUnfold3.C:71
 testUnfold3.C:72
 testUnfold3.C:73
 testUnfold3.C:74
 testUnfold3.C:75
 testUnfold3.C:76
 testUnfold3.C:77
 testUnfold3.C:78
 testUnfold3.C:79
 testUnfold3.C:80
 testUnfold3.C:81
 testUnfold3.C:82
 testUnfold3.C:83
 testUnfold3.C:84
 testUnfold3.C:85
 testUnfold3.C:86
 testUnfold3.C:87
 testUnfold3.C:88
 testUnfold3.C:89
 testUnfold3.C:90
 testUnfold3.C:91
 testUnfold3.C:92
 testUnfold3.C:93
 testUnfold3.C:94
 testUnfold3.C:95
 testUnfold3.C:96
 testUnfold3.C:97
 testUnfold3.C:98
 testUnfold3.C:99
 testUnfold3.C:100
 testUnfold3.C:101
 testUnfold3.C:102
 testUnfold3.C:103
 testUnfold3.C:104
 testUnfold3.C:105
 testUnfold3.C:106
 testUnfold3.C:107
 testUnfold3.C:108
 testUnfold3.C:109
 testUnfold3.C:110
 testUnfold3.C:111
 testUnfold3.C:112
 testUnfold3.C:113
 testUnfold3.C:114
 testUnfold3.C:115
 testUnfold3.C:116
 testUnfold3.C:117
 testUnfold3.C:118
 testUnfold3.C:119
 testUnfold3.C:120
 testUnfold3.C:121
 testUnfold3.C:122
 testUnfold3.C:123
 testUnfold3.C:124
 testUnfold3.C:125
 testUnfold3.C:126
 testUnfold3.C:127
 testUnfold3.C:128
 testUnfold3.C:129
 testUnfold3.C:130
 testUnfold3.C:131
 testUnfold3.C:132
 testUnfold3.C:133
 testUnfold3.C:134
 testUnfold3.C:135
 testUnfold3.C:136
 testUnfold3.C:137
 testUnfold3.C:138
 testUnfold3.C:139
 testUnfold3.C:140
 testUnfold3.C:141
 testUnfold3.C:142
 testUnfold3.C:143
 testUnfold3.C:144
 testUnfold3.C:145
 testUnfold3.C:146
 testUnfold3.C:147
 testUnfold3.C:148
 testUnfold3.C:149
 testUnfold3.C:150
 testUnfold3.C:151
 testUnfold3.C:152
 testUnfold3.C:153
 testUnfold3.C:154
 testUnfold3.C:155
 testUnfold3.C:156
 testUnfold3.C:157
 testUnfold3.C:158
 testUnfold3.C:159
 testUnfold3.C:160
 testUnfold3.C:161
 testUnfold3.C:162
 testUnfold3.C:163
 testUnfold3.C:164
 testUnfold3.C:165
 testUnfold3.C:166
 testUnfold3.C:167
 testUnfold3.C:168
 testUnfold3.C:169
 testUnfold3.C:170
 testUnfold3.C:171
 testUnfold3.C:172
 testUnfold3.C:173
 testUnfold3.C:174
 testUnfold3.C:175
 testUnfold3.C:176
 testUnfold3.C:177
 testUnfold3.C:178
 testUnfold3.C:179
 testUnfold3.C:180
 testUnfold3.C:181
 testUnfold3.C:182
 testUnfold3.C:183
 testUnfold3.C:184
 testUnfold3.C:185
 testUnfold3.C:186
 testUnfold3.C:187
 testUnfold3.C:188
 testUnfold3.C:189
 testUnfold3.C:190
 testUnfold3.C:191
 testUnfold3.C:192
 testUnfold3.C:193
 testUnfold3.C:194
 testUnfold3.C:195
 testUnfold3.C:196
 testUnfold3.C:197
 testUnfold3.C:198
 testUnfold3.C:199
 testUnfold3.C:200
 testUnfold3.C:201
 testUnfold3.C:202
 testUnfold3.C:203
 testUnfold3.C:204
 testUnfold3.C:205
 testUnfold3.C:206
 testUnfold3.C:207
 testUnfold3.C:208
 testUnfold3.C:209
 testUnfold3.C:210
 testUnfold3.C:211
 testUnfold3.C:212
 testUnfold3.C:213
 testUnfold3.C:214
 testUnfold3.C:215
 testUnfold3.C:216
 testUnfold3.C:217
 testUnfold3.C:218
 testUnfold3.C:219
 testUnfold3.C:220
 testUnfold3.C:221
 testUnfold3.C:222
 testUnfold3.C:223
 testUnfold3.C:224
 testUnfold3.C:225
 testUnfold3.C:226
 testUnfold3.C:227
 testUnfold3.C:228
 testUnfold3.C:229
 testUnfold3.C:230
 testUnfold3.C:231
 testUnfold3.C:232
 testUnfold3.C:233
 testUnfold3.C:234
 testUnfold3.C:235
 testUnfold3.C:236
 testUnfold3.C:237
 testUnfold3.C:238
 testUnfold3.C:239
 testUnfold3.C:240
 testUnfold3.C:241
 testUnfold3.C:242
 testUnfold3.C:243
 testUnfold3.C:244
 testUnfold3.C:245
 testUnfold3.C:246
 testUnfold3.C:247
 testUnfold3.C:248
 testUnfold3.C:249
 testUnfold3.C:250
 testUnfold3.C:251
 testUnfold3.C:252
 testUnfold3.C:253
 testUnfold3.C:254
 testUnfold3.C:255
 testUnfold3.C:256
 testUnfold3.C:257
 testUnfold3.C:258
 testUnfold3.C:259
 testUnfold3.C:260
 testUnfold3.C:261
 testUnfold3.C:262
 testUnfold3.C:263
 testUnfold3.C:264
 testUnfold3.C:265
 testUnfold3.C:266
 testUnfold3.C:267
 testUnfold3.C:268
 testUnfold3.C:269
 testUnfold3.C:270
 testUnfold3.C:271
 testUnfold3.C:272
 testUnfold3.C:273
 testUnfold3.C:274
 testUnfold3.C:275
 testUnfold3.C:276
 testUnfold3.C:277
 testUnfold3.C:278
 testUnfold3.C:279
 testUnfold3.C:280
 testUnfold3.C:281
 testUnfold3.C:282
 testUnfold3.C:283
 testUnfold3.C:284
 testUnfold3.C:285
 testUnfold3.C:286
 testUnfold3.C:287
 testUnfold3.C:288
 testUnfold3.C:289
 testUnfold3.C:290
 testUnfold3.C:291
 testUnfold3.C:292
 testUnfold3.C:293
 testUnfold3.C:294
 testUnfold3.C:295
 testUnfold3.C:296
 testUnfold3.C:297
 testUnfold3.C:298
 testUnfold3.C:299
 testUnfold3.C:300
 testUnfold3.C:301
 testUnfold3.C:302
 testUnfold3.C:303
 testUnfold3.C:304
 testUnfold3.C:305
 testUnfold3.C:306
 testUnfold3.C:307
 testUnfold3.C:308
 testUnfold3.C:309
 testUnfold3.C:310
 testUnfold3.C:311
 testUnfold3.C:312
 testUnfold3.C:313
 testUnfold3.C:314
 testUnfold3.C:315
 testUnfold3.C:316
 testUnfold3.C:317
 testUnfold3.C:318
 testUnfold3.C:319
 testUnfold3.C:320
 testUnfold3.C:321
 testUnfold3.C:322
 testUnfold3.C:323
 testUnfold3.C:324
 testUnfold3.C:325
 testUnfold3.C:326
 testUnfold3.C:327
 testUnfold3.C:328
 testUnfold3.C:329
 testUnfold3.C:330
 testUnfold3.C:331
 testUnfold3.C:332
 testUnfold3.C:333
 testUnfold3.C:334
 testUnfold3.C:335
 testUnfold3.C:336
 testUnfold3.C:337
 testUnfold3.C:338
 testUnfold3.C:339
 testUnfold3.C:340
 testUnfold3.C:341
 testUnfold3.C:342
 testUnfold3.C:343
 testUnfold3.C:344
 testUnfold3.C:345
 testUnfold3.C:346
 testUnfold3.C:347
 testUnfold3.C:348
 testUnfold3.C:349
 testUnfold3.C:350
 testUnfold3.C:351
 testUnfold3.C:352
 testUnfold3.C:353
 testUnfold3.C:354
 testUnfold3.C:355
 testUnfold3.C:356
 testUnfold3.C:357
 testUnfold3.C:358
 testUnfold3.C:359
 testUnfold3.C:360
 testUnfold3.C:361
 testUnfold3.C:362
 testUnfold3.C:363
 testUnfold3.C:364
 testUnfold3.C:365
 testUnfold3.C:366
 testUnfold3.C:367
 testUnfold3.C:368
 testUnfold3.C:369
 testUnfold3.C:370
 testUnfold3.C:371
 testUnfold3.C:372
 testUnfold3.C:373
 testUnfold3.C:374
 testUnfold3.C:375
 testUnfold3.C:376
 testUnfold3.C:377
 testUnfold3.C:378
 testUnfold3.C:379
 testUnfold3.C:380
 testUnfold3.C:381
 testUnfold3.C:382
 testUnfold3.C:383
 testUnfold3.C:384
 testUnfold3.C:385
 testUnfold3.C:386
 testUnfold3.C:387
 testUnfold3.C:388
 testUnfold3.C:389
 testUnfold3.C:390
 testUnfold3.C:391
 testUnfold3.C:392
 testUnfold3.C:393
 testUnfold3.C:394
 testUnfold3.C:395
 testUnfold3.C:396
 testUnfold3.C:397
 testUnfold3.C:398
 testUnfold3.C:399
 testUnfold3.C:400
 testUnfold3.C:401
 testUnfold3.C:402
 testUnfold3.C:403
 testUnfold3.C:404
 testUnfold3.C:405
 testUnfold3.C:406
 testUnfold3.C:407
 testUnfold3.C:408
 testUnfold3.C:409
 testUnfold3.C:410
 testUnfold3.C:411
 testUnfold3.C:412
 testUnfold3.C:413
 testUnfold3.C:414
 testUnfold3.C:415
 testUnfold3.C:416
 testUnfold3.C:417
 testUnfold3.C:418
 testUnfold3.C:419
 testUnfold3.C:420
 testUnfold3.C:421
 testUnfold3.C:422
 testUnfold3.C:423
 testUnfold3.C:424
 testUnfold3.C:425
 testUnfold3.C:426
 testUnfold3.C:427
 testUnfold3.C:428
 testUnfold3.C:429
 testUnfold3.C:430
 testUnfold3.C:431
 testUnfold3.C:432
 testUnfold3.C:433
 testUnfold3.C:434
 testUnfold3.C:435
 testUnfold3.C:436
 testUnfold3.C:437
 testUnfold3.C:438
 testUnfold3.C:439
 testUnfold3.C:440
 testUnfold3.C:441
 testUnfold3.C:442
 testUnfold3.C:443
 testUnfold3.C:444
 testUnfold3.C:445
 testUnfold3.C:446
 testUnfold3.C:447
 testUnfold3.C:448
 testUnfold3.C:449
 testUnfold3.C:450
 testUnfold3.C:451
 testUnfold3.C:452
 testUnfold3.C:453
 testUnfold3.C:454
 testUnfold3.C:455
 testUnfold3.C:456
 testUnfold3.C:457
 testUnfold3.C:458
 testUnfold3.C:459
 testUnfold3.C:460
 testUnfold3.C:461
 testUnfold3.C:462
 testUnfold3.C:463
 testUnfold3.C:464
 testUnfold3.C:465
 testUnfold3.C:466
 testUnfold3.C:467
 testUnfold3.C:468
 testUnfold3.C:469
 testUnfold3.C:470
 testUnfold3.C:471
 testUnfold3.C:472
 testUnfold3.C:473
 testUnfold3.C:474
 testUnfold3.C:475
 testUnfold3.C:476
 testUnfold3.C:477
 testUnfold3.C:478
 testUnfold3.C:479
 testUnfold3.C:480
 testUnfold3.C:481
 testUnfold3.C:482
 testUnfold3.C:483
 testUnfold3.C:484
 testUnfold3.C:485
 testUnfold3.C:486
 testUnfold3.C:487
 testUnfold3.C:488
 testUnfold3.C:489
 testUnfold3.C:490
 testUnfold3.C:491
 testUnfold3.C:492
 testUnfold3.C:493
 testUnfold3.C:494
 testUnfold3.C:495
 testUnfold3.C:496
 testUnfold3.C:497
 testUnfold3.C:498
 testUnfold3.C:499
 testUnfold3.C:500
 testUnfold3.C:501
 testUnfold3.C:502
 testUnfold3.C:503
 testUnfold3.C:504
 testUnfold3.C:505
 testUnfold3.C:506
 testUnfold3.C:507
 testUnfold3.C:508
 testUnfold3.C:509
 testUnfold3.C:510
 testUnfold3.C:511
 testUnfold3.C:512
 testUnfold3.C:513
 testUnfold3.C:514
 testUnfold3.C:515
 testUnfold3.C:516
 testUnfold3.C:517
 testUnfold3.C:518
 testUnfold3.C:519
 testUnfold3.C:520
 testUnfold3.C:521
 testUnfold3.C:522
 testUnfold3.C:523
 testUnfold3.C:524
 testUnfold3.C:525
 testUnfold3.C:526
 testUnfold3.C:527
 testUnfold3.C:528
 testUnfold3.C:529
 testUnfold3.C:530
 testUnfold3.C:531
 testUnfold3.C:532
 testUnfold3.C:533
 testUnfold3.C:534
 testUnfold3.C:535
 testUnfold3.C:536
 testUnfold3.C:537
 testUnfold3.C:538
 testUnfold3.C:539
 testUnfold3.C:540
 testUnfold3.C:541
 testUnfold3.C:542
 testUnfold3.C:543
 testUnfold3.C:544
 testUnfold3.C:545
 testUnfold3.C:546
 testUnfold3.C:547
 testUnfold3.C:548
 testUnfold3.C:549
 testUnfold3.C:550
 testUnfold3.C:551
 testUnfold3.C:552
 testUnfold3.C:553
 testUnfold3.C:554
 testUnfold3.C:555
 testUnfold3.C:556
 testUnfold3.C:557
 testUnfold3.C:558
 testUnfold3.C:559
 testUnfold3.C:560
 testUnfold3.C:561
 testUnfold3.C:562
 testUnfold3.C:563
 testUnfold3.C:564
 testUnfold3.C:565
 testUnfold3.C:566
 testUnfold3.C:567
 testUnfold3.C:568
 testUnfold3.C:569
 testUnfold3.C:570
 testUnfold3.C:571
 testUnfold3.C:572
 testUnfold3.C:573
 testUnfold3.C:574
 testUnfold3.C:575
 testUnfold3.C:576
 testUnfold3.C:577
 testUnfold3.C:578
 testUnfold3.C:579
 testUnfold3.C:580
 testUnfold3.C:581
 testUnfold3.C:582
 testUnfold3.C:583
 testUnfold3.C:584
 testUnfold3.C:585
 testUnfold3.C:586
 testUnfold3.C:587
 testUnfold3.C:588
 testUnfold3.C:589
 testUnfold3.C:590
 testUnfold3.C:591
 testUnfold3.C:592
 testUnfold3.C:593
 testUnfold3.C:594
 testUnfold3.C:595
 testUnfold3.C:596
 testUnfold3.C:597
 testUnfold3.C:598
 testUnfold3.C:599
 testUnfold3.C:600
 testUnfold3.C:601
 testUnfold3.C:602
 testUnfold3.C:603
 testUnfold3.C:604
 testUnfold3.C:605
 testUnfold3.C:606
 testUnfold3.C:607
 testUnfold3.C:608
 testUnfold3.C:609
 testUnfold3.C:610
 testUnfold3.C:611
 testUnfold3.C:612
 testUnfold3.C:613
 testUnfold3.C:614
 testUnfold3.C:615
 testUnfold3.C:616
 testUnfold3.C:617
 testUnfold3.C:618
 testUnfold3.C:619
 testUnfold3.C:620
 testUnfold3.C:621
 testUnfold3.C:622
 testUnfold3.C:623
 testUnfold3.C:624
 testUnfold3.C:625
 testUnfold3.C:626
 testUnfold3.C:627
 testUnfold3.C:628