testUnfold3.C: Test program for the class TUnfoldSys

// Test program for the class TUnfoldSys
// Author: Stefan Schmitt
// DESY, 14.10.2008

//  Version 16, parallel to changes in TUnfold
//
//  History:
//     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 <TUnfoldSys.h>

using namespace std;

///////////////////////////////////////////////////////////////////////
// 
// Test program for the class TUnfoldSys
//
//  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
//
//  two types of systematic error are studied:
//    SYS1: variation of the input shape
//    SYS2: variation of the trigger efficiency at threshold
//
//  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 "triggered" 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;

  // Binning
  // reconstructed mass
  Int_t const nDet=24;
  Double_t const xminDet=4.0;
  Double_t const xmaxDet=28.0;
  // signal contributions
  Int_t const nGen=10;
  Double_t const xminGen= 6.0;
  Double_t const xmaxGen=26.0;

  //========================
  // Step 1: fill histograms

  // =================== data true 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%
  };

  //============================================
  // generate data distribution
  TH1D *histDetDATA=new TH1D("PT(data)",";Pt(obs)",nDet,xminDet,xmaxDet);
  // second data distribution, with generator-level bins
  // for bin-by-bin correction study
  TH1D *histDetDATAbbb=new TH1D("PT(data,bbb)",";Pt(obs,bbb)",
                                nGen,xminGen,xmaxGen);
  // in real life we do not have this distribution...
  TH1D *histGenDATA=new TH1D("PT(MC,signal,gen)",";Pt(gen)",
                             nGen,xminGen,xmaxGen);

  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) {
        histDetDATA->Fill(ptObs);
        histDetDATAbbb->Fill(ptObs);
     }
     // fill generator-level signal histogram for data
     // for real data, this does not exist
     if(iTypeGen==0) histGenDATA->Fill(ptGen);
  }


  // =================== 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 !!! UNKNOWN !!! 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 os 95%
  };

  
  // study trigger systematics: change parameters for trigger threshold
  static Double_t triggerParm_MCSYS1[]={7.7, // threshold is 7 GeV
                                        3.7, // width is 3 GeV
                                        0.95 // high Pt efficiency is 9%
  };
  // study dependence on initial signal shape
  static Double_t parm_MC_SYS2[]={
     0.0, // fraction of background: not needed
     0.0, // fraction of background: not needed
     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
  };

  //============================================
  // generate MC distributions
  // detector-level histograms
  TH1D *histDetMC[3];
  histDetMC[0]=new TH1D("PT(MC,signal)",";Pt(obs)",nDet,xminDet,xmaxDet);
  histDetMC[1]=new TH1D("PT(MC,bgr1)",";Pt(obs)",nDet,xminDet,xmaxDet);
  histDetMC[2]=new TH1D("PT(MC,bgr2)",";Pt(obs)",nDet,xminDet,xmaxDet);
  TH1D *histDetMCall=new TH1D("PT(MC,all)",";Pt(obs)",nDet,xminDet,xmaxDet);
  TH1D *histDetMCbgr=new TH1D("PT(MC,bgr)",";Pt(obs)",nDet,xminDet,xmaxDet);
  // another set of detector level histograms, this time with
  // generator-level binning, to try the bin-by-bin correction
  TH1D *histDetMCbbb[3];
  histDetMCbbb[0]=new TH1D("PT(MC,sig)",";Pt(obs,bbb)",nGen,xminGen,xmaxGen);
  histDetMCbbb[1]=new TH1D("PT(MC,bgr1)bbb",";Pt(obs,bbb)",nGen,xminGen,xmaxGen);
  histDetMCbbb[2]=new TH1D("PT(MC,bgr2)bbb",";Pt(obs,bbb)",nGen,xminGen,xmaxGen);
  // true signal distribution, two different binnings
  TH1D *histGenMC=new TH1D("PT(MC,signal,gen)bbb",";Pt(gen)",
                         nGen,xminGen,xmaxGen);
  TH2D *histGenDet=new TH2D("PT(MC,signal,gen,det)",";Pt(gen);Pt(det)",
                            nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);

  // 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);
     // detector-level distributions
     // require trigger
     if(isTriggered) {
        // fill signal, bgr1,bgr2 into separate histograms
        histDetMC[iTypeGen]->Fill(ptObs,lumiWeight);
        histDetMCbbb[iTypeGen]->Fill(ptObs,lumiWeight);
        // fill total MC prediction (signal+bgr1+bgr2)
        histDetMCall->Fill(ptObs,lumiWeight);
        // fill bgr only MC prediction (bgr1+bgr2)
        if(iTypeGen>0) {
           histDetMCbgr->Fill(ptObs,lumiWeight);
        }
     }
     // generator level distributions (signal only)
     if(iTypeGen==0) {
        histGenMC->Fill(ptGen,lumiWeight);
     }
     // matrix dectribing the migrations (signal only)
     if(iTypeGen==0) {
        // if event was triggered, fill histogram with ptObs
        if(isTriggered) {
           histGenDet->Fill(ptGen,ptObs,lumiWeight);
        } else {
           // not triggered: fill detector-level underflow bin
           histGenDet->Fill(ptGen,-999.,lumiWeight);
        }
     }
  }

  // generate another MC, this time with changed trigger settings
  // fill 2D histogram and histograms for bin-by-bin study
  TH2D *histGenDetSYS1=new TH2D("PT(MC,signal,gen,det,sys1)",
                                ";Pt(gen);Pt(obs)",
                                nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
  TH1D *histGenSYS1=new TH1D("PT(MC,signal,gen,sys1)",
                             ";Pt(obs)",
                             nGen,xminGen,xmaxGen);
  TH1D *histDetSYS1bbb=new TH1D("PT(MC,signal,det,sys1,bbb)",
                                ";Pt(obs)",
                                nGen,xminGen,xmaxGen);
  intLumi=0.;
  while(intLumi<lumiMC) {
     Int_t iTypeGen;
     Bool_t isTriggered;
     Double_t ptGen;
     Double_t ptObs=GenerateEvent(parm_MC,triggerParm_MCSYS1,&intLumi,
                                  &isTriggered,&ptGen,&iTypeGen);
     // matrix describing the migrations (signal only)
     if(iTypeGen==0) {
        // if event was triggered, fill histogram with ptObs
        if(isTriggered) {
           histGenDetSYS1->Fill(ptGen,ptObs,lumiWeight);
        } else {
           // not triggered: fill detector-level underflow bin
           histGenDetSYS1->Fill(ptGen,-999.,lumiWeight);
        }
     }
     // generator level distribution for bin-by-bin study
     if(iTypeGen==0) {
        histGenSYS1->Fill(ptGen,lumiWeight);
        // detector level for bin-by-bin study
        if(isTriggered) {
           histDetSYS1bbb->Fill(ptObs,lumiWeight);
        }
     }
  }
  // generate a third MC, this time with changed initial distribution
  // only the 2D histogram is filled
  TH2D *histGenDetSYS2=new TH2D("PT(MC,signal,gen,det,sys2)",
                                ";Pt(gen);Pt(det)",
                                nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
  TH1D *histGenSYS2=new TH1D("PT(MC,signal,gen,sys2)",
                             ";Pt(obs)",
                             nGen,xminGen,xmaxGen);
  TH1D *histDetSYS2bbb=new TH1D("PT(MC,signal,det,sys2,bbb)",
                                ";Pt(obs)",
                                nGen,xminGen,xmaxGen);
  intLumi=0.0;
  while(intLumi<lumiMC) {
     Int_t iTypeGen;
     Bool_t isTriggered;
     Double_t ptGen;
     Double_t ptObs=GenerateEvent(parm_MC_SYS2,triggerParm_MC,&intLumi,
                                  &isTriggered,&ptGen,&iTypeGen);
     // matrix dectribing the migrations (signal only)
     if(iTypeGen==0) {
        // if event was triggered, fill histogram with ptObs
        if(isTriggered) {
           histGenDetSYS2->Fill(ptGen,ptObs,lumiWeight);
        } else {
           // not triggered: fill detector-level underflow bin
           histGenDetSYS2->Fill(ptGen,-999.,lumiWeight);
        }
     }
     if(iTypeGen==0) {
        // generator level distribution for bin-by-bin study
        histGenSYS2->Fill(ptGen,lumiWeight);
        // detector level for bin-by-bin study
        if(isTriggered) {
           histDetSYS2bbb->Fill(ptObs,lumiWeight);
        }
     }
  }

  //========================
  // Step 2: unfolding
  //
  // this is based on a matrix (TH2D) which describes the connection
  // between
  //   nDet    Detector bins
  //   nGen+2  Generator level bins
  //
  // why are there nGen+2 Generator level bins?
  //   the two extra bins are the underflow and overflow bin
  //   These are unfolded as well. Such bins are needed to
  //   accomodate for migrations from outside the phase-space into the
  //   observed phase-space.
  //
  // In addition there are overflow/underflow bins
  //   for the reconstructed variable. These are used to count events
  //   which are -NOT- reconstructed. Either because they migrate
  //   out of the reconstructed phase-space. Or because there are detector
  //   inefficiencies.
  //
  TUnfoldSys unfold(histGenDet,TUnfold::kHistMapOutputHoriz,
                    TUnfold::kRegModeSize,TUnfold::kEConstraintNone);

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

  // subtract background, normalized to data luminosity
  //  with 10% scale error each
  Double_t scale_bgr=1.0;
  Double_t dscale_bgr=0.2;
  unfold.SubtractBackground(histDetMC[1],"background1",scale_bgr,dscale_bgr);
  unfold.SubtractBackground(histDetMC[2],"background2",scale_bgr,dscale_bgr);

  // add systematic errors
  // trigger threshold
  unfold.AddSysError(histGenDetSYS1,"trigger_SYS1",
                     TUnfold::kHistMapOutputHoriz,
                     TUnfoldSys::kSysErrModeMatrix);
  // calorimeter response
  unfold.AddSysError(histGenDetSYS2,"signalshape_SYS2",
                     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);

  // 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);

  // get unfolding output
  // Reminder: there are
  //   nGen+2  Generator level bins
  // But we are not interested in the additional bins.
  // The mapping is required to get the proper errors and correlations

  // Here a mapping is defined from the nGen+2 bins to nGen bins.
  Int_t *binMap=new Int_t[nGen+2];
  binMap[0] = -1;    // underflow bin should be ignored
  binMap[nGen+1]=-1; // overflow bin should be ignored
  for(Int_t i=1;i<=nGen;i++) {
     binMap[i]=i;  // map inner bins to output histogram bins
  }

  // this returns the unfolding output
  // The errors included at this point are:
  //    * statistical errros
  //    * background subtraction errors
  TH1D *histUnfoldStatBgr=new TH1D("PT(unfold,signal,stat+bgr)",";Pt(gen)",
                            nGen,xminGen,xmaxGen);
  unfold.GetOutput(histUnfoldStatBgr,binMap);

  // retreive error matrix of statistical errors only
  TH2D *histEmatStat=new TH2D("ErrMat(stat)",";Pt(gen);Pt(gen)",
                              nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
  unfold.GetEmatrixInput(histEmatStat,binMap);

  // retreive full error matrix
  // This includes also all systematic errors defined above
  TH2D *histEmatTotal=new TH2D("ErrMat(total)",";Pt(gen);Pt(gen)",
                               nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
  unfold.GetEmatrixTotal(histEmatTotal,binMap);

  // create two copies of the unfolded data, one with statistical errors
  // the other with total errors
  TH1D *histUnfoldStat=new TH1D("PT(unfold,signal,stat)",";Pt(gen)",
                                nGen,xminGen,xmaxGen);
  TH1D *histUnfoldTotal=new TH1D("PT(unfold,signal,total)",";Pt(gen)",
                                 nGen,xminGen,xmaxGen);
  for(Int_t i=0;i<nGen+2;i++) {
     Double_t c=histUnfoldStatBgr->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);
     }
  }

  delete [] binMap;

  //========================
  // Step 3: plots

  TCanvas output;
  output.Divide(3,2);
  output.cd(1);
  histDetDATA->SetMinimum(0.0);
  histDetDATA->Draw("E");
  histDetMCall->SetMinimum(0.0);
  histDetMCall->SetLineColor(kBlue);  
  histDetMCbgr->SetLineColor(kRed);  
  histDetMC[1]->SetLineColor(kCyan);  
  histDetMCall->Draw("SAME HIST");
  histDetMCbgr->Draw("SAME HIST");
  histDetMC[1]->Draw("SAME HIST");

  output.cd(2);
  histUnfoldTotal->SetMinimum(0.0);
  histUnfoldTotal->SetMaximum(histUnfoldTotal->GetMaximum()*1.5);
  // outer error: total error
  histUnfoldTotal->Draw("E");
  // middle error: stat+bgr
  histUnfoldStatBgr->Draw("SAME E1");
  // inner error: stat only
  histUnfoldStat->Draw("SAME E1");
  histGenDATA->Draw("SAME HIST");
  histGenMC->SetLineColor(kBlue);
  histGenMC->Draw("SAME HIST");

  output.cd(3);
  histGenDet->SetLineColor(kBlue);
  histGenDet->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 4: compare results to
  // the so-called bin-by-bin "correction"
  for(Int_t i=1;i<=nGen;i++) {
     // data contribution in this bin
     Double_t data=histDetDATAbbb->GetBinContent(i);
     Double_t errData=histDetDATAbbb->GetBinError(i);

     // subtract background contribution
     Double_t data_bgr=data;
     Double_t errData_bgr=errData;
     for(Int_t j=1;j<=2;j++) {
        data_bgr -= scale_bgr*histDetMCbbb[j]->GetBinContent(i);
        errData_bgr = TMath::Sqrt(errData_bgr*errData_bgr+
                                  scale_bgr*histDetMCbbb[j]->GetBinError(i)*
                                  scale_bgr*histDetMCbbb[j]->GetBinError(i)+
                                  dscale_bgr*histDetMCbbb[j]->GetBinContent(i)*
                                  dscale_bgr*histDetMCbbb[j]->GetBinContent(i)
                                  );
     }
     // "correct" the data, using the Monte Carlo and neglecting off-diagonals
     Double_t fCorr=(histGenMC->GetBinContent(i)/
                      histDetMCbbb[0]->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_1=(histGenSYS1->GetBinContent(i)/
                       histDetSYS1bbb->GetBinContent(i));
     Double_t shift_sys1= data_bgr*fCorr_1 - data_bbb;
     Double_t fCorr_2=(histGenSYS2->GetBinContent(i)/
                        histDetSYS2bbb->GetBinContent(i));
     Double_t shift_sys2= data_bgr*fCorr_2 - data_bbb;

     cout<<data_bbb<<" "<<shift_sys1<<" "<<shift_sys2<<"\n";

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

     // 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=histUnfoldStatBgr->GetBinError(i);
     Double_t errData_total_unfold=histUnfoldTotal->GetBinError(i);

     // compare
     std::cout<<"Bin "<<i<<": true "<<histGenDATA->GetBinContent(i)
              <<" unfold: "<<data_unfold
              <<" +/- "<<errData_stat_unfold<<" (stat)"
              <<" +/- "<<TMath::Sqrt(errData_statbgr_unfold*
                                     errData_statbgr_unfold-
                                     errData_stat_unfold*
                                     errData_stat_unfold)<<" (bgr)"
              <<" +/- "<<TMath::Sqrt(errData_total_unfold*
                                     errData_total_unfold-
                                     errData_statbgr_unfold*
                                     errData_statbgr_unfold)<<" (sys)"<<"\n";
     std::cout<<"Bin "<<i<<": true "<<histGenDATA->GetBinContent(i)
              <<" binbybin: "<<data_bbb
              <<" +/- "<<errData_stat_bbb<<" (stat)"
              <<" +/- "<<TMath::Sqrt(errData_statbgr_bbb*
                                     errData_statbgr_bbb-
                                     errData_stat_bbb*
                                     errData_stat_bbb)<<" (bgr)"
              <<" +/- "<<TMath::Sqrt(errData_total_bbb*
                                     errData_total_bbb-
                                     errData_statbgr_bbb*
                                     errData_statbgr_bbb)<<" (sys)"
              <<"\n";
  }
}