From $ROOTSYS/tutorials/math/testUnfold2.C

// Author: Stefan Schmitt
// DESY, 14.10.2008

//  Version 17.0, updated for changed methods in TUnfold
//
//  History:
//    Version 16.1, parallel to changes in TUnfold
//    Version 16.0, parallel to changes in TUnfold
//    Version 15, with automatic L-curve scan, simplified example
//    Version 14, with changes in TUnfoldSys.cxx
//    Version 13,  with changes to TUnfold.C
//    Version 12,  with improvements to TUnfold.cxx
//    Version 11,  print chi**2 and number of degrees of freedom
//    Version 10, with bug-fix in TUnfold.cxx
//    Version 9, with bug-fix in TUnfold.cxx, TUnfold.h
//    Version 8, with bug-fix in TUnfold.cxx, TUnfold.h
//    Version 7, with bug-fix in TUnfold.cxx, TUnfold.h
//    Version 6a, fix problem with dynamic array allocation under windows
//    Version 6, re-include class MyUnfold in the example
//    Version 5, move class MyUnfold to seperate files
//    Version 4, with bug-fix in TUnfold.C
//    Version 3, with bug-fix in TUnfold.C
//    Version 2, with changed ScanLcurve() arguments
//    Version 1, remove L curve analysis, use ScanLcurve() method instead
//    Version 0, L curve analysis included here

/*
  This file is part of TUnfold.

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

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

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

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

using namespace std;

///////////////////////////////////////////////////////////////////////
//
//  Test program as an example for a user specific regularisation scheme
//
//  (1) Generate Monte Carlo and Data events
//      The events consist of
//        signal
//        background
//
//      The signal is a resonance. It is generated with a Breit-Wigner,
//      smeared by a Gaussian
//
//  (2) Unfold the data. The result is:
//      The background level
//      The shape of the resonance, corrected for detector effects
//
//      The regularisation is done on the curvature, excluding the bins
//      near the peak.
//
//  (3) produce some plots
//
///////////////////////////////////////////////////////////////////////

TRandom *rnd=0;

// generate an event
// output:
//  negative mass: background event
//  positive mass: signal event
Double_t GenerateEvent(Double_t bgr, // relative fraction of background
                       Double_t mass, // peak position
                       Double_t gamma) // peak width
{
  Double_t t;
  if(rnd->Rndm()>bgr) {
    // generate signal event
    // with positive mass
    do {
      do {
        t=rnd->Rndm();
      } while(t>=1.0);
      t=TMath::Tan((t-0.5)*TMath::Pi())*gamma+mass;
    } while(t<=0.0);
    return t;
  } else {
    // generate background event
    // generate events following a power-law distribution
    //   f(E) = K * TMath::power((E0+E),N0)
    static Double_t const E0=2.4;
    static Double_t const N0=2.9;
    do {
      do {
        t=rnd->Rndm();
      } while(t>=1.0);
      // the mass is returned negative
      // In our example a convenient way to indicate it is a background event.
      t= -(TMath::Power(1.-t,1./(1.-N0))-1.0)*E0;
    } while(t>=0.0);
    return t;
  }
}

// smear the event to detector level
// input:
//   mass on generator level (mTrue>0 !)
// output:
//   mass on detector level
Double_t DetectorEvent(Double_t mTrue) {
  // smear by double-gaussian
  static Double_t frac=0.1;
  static Double_t wideBias=0.03;
  static Double_t wideSigma=0.5;
  static Double_t smallBias=0.0;
  static Double_t smallSigma=0.1;
  if(rnd->Rndm()>frac) {
    return rnd->Gaus(mTrue+smallBias,smallSigma);
  } else {
    return rnd->Gaus(mTrue+wideBias,wideSigma);
  }
}

int testUnfold2()
{
  // switch on histogram errors
  TH1::SetDefaultSumw2();

  // random generator
  rnd=new TRandom3();

  // data and MC luminosity, cross-section
  Double_t const luminosityData=100000;
  Double_t const luminosityMC=1000000;
  Double_t const crossSection=1.0;

  Int_t const nDet=250;
  Int_t const nGen=100;
  Double_t const xminDet=0.0;
  Double_t const xmaxDet=10.0;
  Double_t const xminGen=0.0;
  Double_t const xmaxGen=10.0;

  //============================================
  // generate MC distribution
  //
  TH1D *histMgenMC=new TH1D("MgenMC",";mass(gen)",nGen,xminGen,xmaxGen);
  TH1D *histMdetMC=new TH1D("MdetMC",";mass(det)",nDet,xminDet,xmaxDet);
  TH2D *histMdetGenMC=new TH2D("MdetgenMC",";mass(det);mass(gen)",nDet,xminDet,xmaxDet,
                              nGen,xminGen,xmaxGen);
  Int_t neventMC=rnd->Poisson(luminosityMC*crossSection);
  for(Int_t i=0;i<neventMC;i++) {
    Double_t mGen=GenerateEvent(0.3, // relative fraction of background
                                4.0, // peak position in MC
                                0.2); // peak width in MC
    Double_t mDet=DetectorEvent(TMath::Abs(mGen));
    // the generated mass is negative for background
    // and positive for signal
    // so it will be filled in the underflow bin
    // this is very convenient for the unfolding:
    // the unfolded result will contain the number of background
    // events in the underflow bin

    // generated MC distribution (for comparison only)
    histMgenMC->Fill(mGen,luminosityData/luminosityMC);
    // reconstructed MC distribution (for comparison only)
    histMdetMC->Fill(mDet,luminosityData/luminosityMC);

    // matrix describing how the generator input migrates to the
    // reconstructed level. Unfolding input.
    // NOTE on underflow/overflow bins:
    //  (1) the detector level under/overflow bins are used for
    //       normalisation ("efficiency" correction)
    //       in our toy example, these bins are populated from tails
    //       of the initial MC distribution.
    //  (2) the generator level underflow/overflow bins are
    //       unfolded. In this example:
    //       underflow bin: background events reconstructed in the detector
    //       overflow bin: signal events generated at masses > xmaxDet
    // for the unfolded result these bins will be filled
    //  -> the background normalisation will be contained in the underflow bin
    histMdetGenMC->Fill(mDet,mGen,luminosityData/luminosityMC);
  }

  //============================================
  // generate data distribution
  //
  TH1D *histMgenData=new TH1D("MgenData",";mass(gen)",nGen,xminGen,xmaxGen);
  TH1D *histMdetData=new TH1D("MdetData",";mass(det)",nDet,xminDet,xmaxDet);
  Int_t neventData=rnd->Poisson(luminosityData*crossSection);
  for(Int_t i=0;i<neventData;i++) {
    Double_t mGen=GenerateEvent(0.4, // relative fraction of background
                                3.8, // peak position
                                0.15); // peak width
    Double_t mDet=DetectorEvent(TMath::Abs(mGen));
    // generated data mass for comparison plots
    // for real data, we do not have this histogram
    histMgenData->Fill(mGen);

    // reconstructed mass, unfolding input
    histMdetData->Fill(mDet);
  }

  //=========================================================================
  // set up the unfolding
  TUnfold unfold(histMdetGenMC,TUnfold::kHistMapOutputVert,
                 TUnfold::kRegModeNone);
  // regularisation
  //----------------
  // the regularisation is done on the curvature (2nd derivative) of
  // the output distribution
  //
  // One has to exclude the bins near the peak of the Breit-Wigner,
  // because there the curvature is high
  // (and the regularisation eventually could enforce a small
  //  curvature, thus biasing result)
  //
  // in real life, the parameters below would have to be optimized,
  // depending on the data peak position and width
  // Or maybe one finds a different regularisation scheme... this is
  // just an example...
  Double_t estimatedPeakPosition=3.8;
  Int_t nPeek=3;
  TUnfold::ERegMode regMode=TUnfold::kRegModeCurvature;
  // calculate bin number correspoinding to estimated peak position
  Int_t iPeek=(Int_t)(nGen*(estimatedPeakPosition-xminGen)/(xmaxGen-xminGen)
                      // offset 1.5
                      // accounts for start bin 1
                      // and rounding errors +0.5
                      +1.5);
  // regularize output bins 1..iPeek-nPeek
  unfold.RegularizeBins(1,1,iPeek-nPeek,regMode);
  // regularize output bins iPeek+nPeek..nGen
  unfold.RegularizeBins(iPeek+nPeek,1,nGen-(iPeek+nPeek),regMode);

  // unfolding
  //-----------

  // set input distribution and bias scale (=0)
  if(unfold.SetInput(histMdetData,0.0)>=10000) {
    std::cout<<"Unfolding result may be wrong\n";
  }

  // do the unfolding here
  Double_t tauMin=0.0;
  Double_t tauMax=0.0;
  Int_t nScan=30;
  Int_t iBest;
  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
  iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
  std::cout<<"tau="<<unfold.GetTau()<<"\n";
  std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
           <<" / "<<unfold.GetNdf()<<"\n";

  // save point corresponding to the kink in the L curve as TGraph
  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 *bestLogTauX=new TGraph(1,t,x);

  //============================================================
  // extract unfolding results into histograms

  // set up a bin map, excluding underflow and overflow bins
  // the binMap relates the the output of the unfolding to the final
  // histogram bins
  Int_t *binMap=new Int_t[nGen+2];
  for(Int_t i=1;i<=nGen;i++) binMap[i]=i;
  binMap[0]=-1;
  binMap[nGen+1]=-1;

  TH1D *histMunfold=new TH1D("Unfolded",";mass(gen)",nGen,xminGen,xmaxGen);
  unfold.GetOutput(histMunfold,binMap);
  TH1D *histMdetFold=new TH1D("FoldedBack","mass(det)",nDet,xminDet,xmaxDet);
  unfold.GetFoldedOutput(histMdetFold);

  // store global correlation coefficients
  TH1D *histRhoi=new TH1D("rho_I","mass",nGen,xminGen,xmaxGen);
  unfold.GetRhoI(histRhoi,binMap);

  delete[] binMap;
  binMap=0;

  //=====================================================================
  // plot some histograms
  TCanvas output;

  // produce some plots
  output.Divide(3,2);

  // Show the matrix which connects input and output
  // There are overflow bins at the bottom, not shown in the plot
  // These contain the background shape.
  // The overflow bins to the left and right contain
  // events which are not reconstructed. These are necessary for proper MC
  // normalisation
  output.cd(1);
  histMdetGenMC->Draw("BOX");

  // draw generator-level distribution:
  //   data (red) [for real data this is not available]
  //   MC input (black) [with completely wrong peak position and shape]
  //   unfolded data (blue)
  output.cd(2);
  histMunfold->SetLineColor(kBlue);
  histMunfold->Draw();
  histMgenData->SetLineColor(kRed);
  histMgenData->Draw("SAME");
  histMgenMC->Draw("SAME HIST");

  // show detector level distributions
  //    data (red)
  //    MC (black)
  //    unfolded data (blue)
  output.cd(3);
  histMdetFold->SetLineColor(kBlue);
  histMdetFold->Draw();
  histMdetData->SetLineColor(kRed);
  histMdetData->Draw("SAME");
  histMdetMC->Draw("SAME HIST");

  // show correlation coefficients
  //     all bins outside the peak are found to be highly correlated
  //     But they are compatible with zero anyway
  //     If the peak shape is fitted,
  //     these correlations have to be taken into account, see example
  output.cd(4);
  histRhoi->Draw();

  // show rhoi_max(tau) distribution
  output.cd(5);
  logTauX->Draw();
  bestLogTauX->SetMarkerColor(kRed);
  bestLogTauX->Draw("*");

  output.cd(6);
  lCurve->Draw("AL");
  bestLcurve->SetMarkerColor(kRed);
  bestLcurve->Draw("*");

  output.SaveAs("testUnfold2.ps");
  return 0;
}
 testUnfold2.C:1
 testUnfold2.C:2
 testUnfold2.C:3
 testUnfold2.C:4
 testUnfold2.C:5
 testUnfold2.C:6
 testUnfold2.C:7
 testUnfold2.C:8
 testUnfold2.C:9
 testUnfold2.C:10
 testUnfold2.C:11
 testUnfold2.C:12
 testUnfold2.C:13
 testUnfold2.C:14
 testUnfold2.C:15
 testUnfold2.C:16
 testUnfold2.C:17
 testUnfold2.C:18
 testUnfold2.C:19
 testUnfold2.C:20
 testUnfold2.C:21
 testUnfold2.C:22
 testUnfold2.C:23
 testUnfold2.C:24
 testUnfold2.C:25
 testUnfold2.C:26
 testUnfold2.C:27
 testUnfold2.C:28
 testUnfold2.C:29
 testUnfold2.C:30
 testUnfold2.C:31
 testUnfold2.C:32
 testUnfold2.C:33
 testUnfold2.C:34
 testUnfold2.C:35
 testUnfold2.C:36
 testUnfold2.C:37
 testUnfold2.C:38
 testUnfold2.C:39
 testUnfold2.C:40
 testUnfold2.C:41
 testUnfold2.C:42
 testUnfold2.C:43
 testUnfold2.C:44
 testUnfold2.C:45
 testUnfold2.C:46
 testUnfold2.C:47
 testUnfold2.C:48
 testUnfold2.C:49
 testUnfold2.C:50
 testUnfold2.C:51
 testUnfold2.C:52
 testUnfold2.C:53
 testUnfold2.C:54
 testUnfold2.C:55
 testUnfold2.C:56
 testUnfold2.C:57
 testUnfold2.C:58
 testUnfold2.C:59
 testUnfold2.C:60
 testUnfold2.C:61
 testUnfold2.C:62
 testUnfold2.C:63
 testUnfold2.C:64
 testUnfold2.C:65
 testUnfold2.C:66
 testUnfold2.C:67
 testUnfold2.C:68
 testUnfold2.C:69
 testUnfold2.C:70
 testUnfold2.C:71
 testUnfold2.C:72
 testUnfold2.C:73
 testUnfold2.C:74
 testUnfold2.C:75
 testUnfold2.C:76
 testUnfold2.C:77
 testUnfold2.C:78
 testUnfold2.C:79
 testUnfold2.C:80
 testUnfold2.C:81
 testUnfold2.C:82
 testUnfold2.C:83
 testUnfold2.C:84
 testUnfold2.C:85
 testUnfold2.C:86
 testUnfold2.C:87
 testUnfold2.C:88
 testUnfold2.C:89
 testUnfold2.C:90
 testUnfold2.C:91
 testUnfold2.C:92
 testUnfold2.C:93
 testUnfold2.C:94
 testUnfold2.C:95
 testUnfold2.C:96
 testUnfold2.C:97
 testUnfold2.C:98
 testUnfold2.C:99
 testUnfold2.C:100
 testUnfold2.C:101
 testUnfold2.C:102
 testUnfold2.C:103
 testUnfold2.C:104
 testUnfold2.C:105
 testUnfold2.C:106
 testUnfold2.C:107
 testUnfold2.C:108
 testUnfold2.C:109
 testUnfold2.C:110
 testUnfold2.C:111
 testUnfold2.C:112
 testUnfold2.C:113
 testUnfold2.C:114
 testUnfold2.C:115
 testUnfold2.C:116
 testUnfold2.C:117
 testUnfold2.C:118
 testUnfold2.C:119
 testUnfold2.C:120
 testUnfold2.C:121
 testUnfold2.C:122
 testUnfold2.C:123
 testUnfold2.C:124
 testUnfold2.C:125
 testUnfold2.C:126
 testUnfold2.C:127
 testUnfold2.C:128
 testUnfold2.C:129
 testUnfold2.C:130
 testUnfold2.C:131
 testUnfold2.C:132
 testUnfold2.C:133
 testUnfold2.C:134
 testUnfold2.C:135
 testUnfold2.C:136
 testUnfold2.C:137
 testUnfold2.C:138
 testUnfold2.C:139
 testUnfold2.C:140
 testUnfold2.C:141
 testUnfold2.C:142
 testUnfold2.C:143
 testUnfold2.C:144
 testUnfold2.C:145
 testUnfold2.C:146
 testUnfold2.C:147
 testUnfold2.C:148
 testUnfold2.C:149
 testUnfold2.C:150
 testUnfold2.C:151
 testUnfold2.C:152
 testUnfold2.C:153
 testUnfold2.C:154
 testUnfold2.C:155
 testUnfold2.C:156
 testUnfold2.C:157
 testUnfold2.C:158
 testUnfold2.C:159
 testUnfold2.C:160
 testUnfold2.C:161
 testUnfold2.C:162
 testUnfold2.C:163
 testUnfold2.C:164
 testUnfold2.C:165
 testUnfold2.C:166
 testUnfold2.C:167
 testUnfold2.C:168
 testUnfold2.C:169
 testUnfold2.C:170
 testUnfold2.C:171
 testUnfold2.C:172
 testUnfold2.C:173
 testUnfold2.C:174
 testUnfold2.C:175
 testUnfold2.C:176
 testUnfold2.C:177
 testUnfold2.C:178
 testUnfold2.C:179
 testUnfold2.C:180
 testUnfold2.C:181
 testUnfold2.C:182
 testUnfold2.C:183
 testUnfold2.C:184
 testUnfold2.C:185
 testUnfold2.C:186
 testUnfold2.C:187
 testUnfold2.C:188
 testUnfold2.C:189
 testUnfold2.C:190
 testUnfold2.C:191
 testUnfold2.C:192
 testUnfold2.C:193
 testUnfold2.C:194
 testUnfold2.C:195
 testUnfold2.C:196
 testUnfold2.C:197
 testUnfold2.C:198
 testUnfold2.C:199
 testUnfold2.C:200
 testUnfold2.C:201
 testUnfold2.C:202
 testUnfold2.C:203
 testUnfold2.C:204
 testUnfold2.C:205
 testUnfold2.C:206
 testUnfold2.C:207
 testUnfold2.C:208
 testUnfold2.C:209
 testUnfold2.C:210
 testUnfold2.C:211
 testUnfold2.C:212
 testUnfold2.C:213
 testUnfold2.C:214
 testUnfold2.C:215
 testUnfold2.C:216
 testUnfold2.C:217
 testUnfold2.C:218
 testUnfold2.C:219
 testUnfold2.C:220
 testUnfold2.C:221
 testUnfold2.C:222
 testUnfold2.C:223
 testUnfold2.C:224
 testUnfold2.C:225
 testUnfold2.C:226
 testUnfold2.C:227
 testUnfold2.C:228
 testUnfold2.C:229
 testUnfold2.C:230
 testUnfold2.C:231
 testUnfold2.C:232
 testUnfold2.C:233
 testUnfold2.C:234
 testUnfold2.C:235
 testUnfold2.C:236
 testUnfold2.C:237
 testUnfold2.C:238
 testUnfold2.C:239
 testUnfold2.C:240
 testUnfold2.C:241
 testUnfold2.C:242
 testUnfold2.C:243
 testUnfold2.C:244
 testUnfold2.C:245
 testUnfold2.C:246
 testUnfold2.C:247
 testUnfold2.C:248
 testUnfold2.C:249
 testUnfold2.C:250
 testUnfold2.C:251
 testUnfold2.C:252
 testUnfold2.C:253
 testUnfold2.C:254
 testUnfold2.C:255
 testUnfold2.C:256
 testUnfold2.C:257
 testUnfold2.C:258
 testUnfold2.C:259
 testUnfold2.C:260
 testUnfold2.C:261
 testUnfold2.C:262
 testUnfold2.C:263
 testUnfold2.C:264
 testUnfold2.C:265
 testUnfold2.C:266
 testUnfold2.C:267
 testUnfold2.C:268
 testUnfold2.C:269
 testUnfold2.C:270
 testUnfold2.C:271
 testUnfold2.C:272
 testUnfold2.C:273
 testUnfold2.C:274
 testUnfold2.C:275
 testUnfold2.C:276
 testUnfold2.C:277
 testUnfold2.C:278
 testUnfold2.C:279
 testUnfold2.C:280
 testUnfold2.C:281
 testUnfold2.C:282
 testUnfold2.C:283
 testUnfold2.C:284
 testUnfold2.C:285
 testUnfold2.C:286
 testUnfold2.C:287
 testUnfold2.C:288
 testUnfold2.C:289
 testUnfold2.C:290
 testUnfold2.C:291
 testUnfold2.C:292
 testUnfold2.C:293
 testUnfold2.C:294
 testUnfold2.C:295
 testUnfold2.C:296
 testUnfold2.C:297
 testUnfold2.C:298
 testUnfold2.C:299
 testUnfold2.C:300
 testUnfold2.C:301
 testUnfold2.C:302
 testUnfold2.C:303
 testUnfold2.C:304
 testUnfold2.C:305
 testUnfold2.C:306
 testUnfold2.C:307
 testUnfold2.C:308
 testUnfold2.C:309
 testUnfold2.C:310
 testUnfold2.C:311
 testUnfold2.C:312
 testUnfold2.C:313
 testUnfold2.C:314
 testUnfold2.C:315
 testUnfold2.C:316
 testUnfold2.C:317
 testUnfold2.C:318
 testUnfold2.C:319
 testUnfold2.C:320
 testUnfold2.C:321
 testUnfold2.C:322
 testUnfold2.C:323
 testUnfold2.C:324
 testUnfold2.C:325
 testUnfold2.C:326
 testUnfold2.C:327
 testUnfold2.C:328
 testUnfold2.C:329
 testUnfold2.C:330
 testUnfold2.C:331
 testUnfold2.C:332
 testUnfold2.C:333
 testUnfold2.C:334
 testUnfold2.C:335
 testUnfold2.C:336
 testUnfold2.C:337
 testUnfold2.C:338
 testUnfold2.C:339
 testUnfold2.C:340
 testUnfold2.C:341
 testUnfold2.C:342
 testUnfold2.C:343
 testUnfold2.C:344
 testUnfold2.C:345
 testUnfold2.C:346
 testUnfold2.C:347
 testUnfold2.C:348
 testUnfold2.C:349
 testUnfold2.C:350
 testUnfold2.C:351
 testUnfold2.C:352
 testUnfold2.C:353
 testUnfold2.C:354
 testUnfold2.C:355
 testUnfold2.C:356
 testUnfold2.C:357
 testUnfold2.C:358
 testUnfold2.C:359
 testUnfold2.C:360
 testUnfold2.C:361