From e43b090418b0729ab174c905a1b9929c499a58da Mon Sep 17 00:00:00 2001
From: Olivier Couet <olivier.couet@cern.ch>
Date: Wed, 23 Nov 2016 15:47:29 +0100
Subject: [PATCH] Add new tutorials for TUnfold. Stefan Schmitt provided the
examples. More work was needed to: - Create the index.md file - Make the
source files are Doxygen Compliant (formatting, Math formulae) - Spell check
---
tutorials/unfold/index.md | 3 +
tutorials/unfold/testUnfold1.C | 517 +++++++
tutorials/unfold/testUnfold2.C | 366 +++++
tutorials/unfold/testUnfold3.C | 631 +++++++++
tutorials/unfold/testUnfold4.C | 222 +++
tutorials/unfold/testUnfold5a.C | 338 +++++
tutorials/unfold/testUnfold5b.C | 182 +++
tutorials/unfold/testUnfold5c.C | 290 ++++
tutorials/unfold/testUnfold5d.C | 296 ++++
tutorials/unfold/testUnfold6.C | 117 ++
tutorials/unfold/testUnfold6binning.xml | 27 +
tutorials/unfold/testUnfold7a.C | 447 ++++++
tutorials/unfold/testUnfold7b.C | 302 +++++
tutorials/unfold/testUnfold7binning.xml | 18 +
tutorials/unfold/testUnfold7c.C | 1655 +++++++++++++++++++++++
15 files changed, 5411 insertions(+)
create mode 100644 tutorials/unfold/index.md
create mode 100644 tutorials/unfold/testUnfold1.C
create mode 100644 tutorials/unfold/testUnfold2.C
create mode 100644 tutorials/unfold/testUnfold3.C
create mode 100644 tutorials/unfold/testUnfold4.C
create mode 100644 tutorials/unfold/testUnfold5a.C
create mode 100644 tutorials/unfold/testUnfold5b.C
create mode 100644 tutorials/unfold/testUnfold5c.C
create mode 100644 tutorials/unfold/testUnfold5d.C
create mode 100644 tutorials/unfold/testUnfold6.C
create mode 100644 tutorials/unfold/testUnfold6binning.xml
create mode 100644 tutorials/unfold/testUnfold7a.C
create mode 100644 tutorials/unfold/testUnfold7b.C
create mode 100644 tutorials/unfold/testUnfold7binning.xml
create mode 100644 tutorials/unfold/testUnfold7c.C
diff --git a/tutorials/unfold/index.md b/tutorials/unfold/index.md
new file mode 100644
index 00000000000..a6a554a2a2a
--- /dev/null
+++ b/tutorials/unfold/index.md
@@ -0,0 +1,3 @@
+\defgroup tutorial_unfold Unfold tutorials
+\ingroup Tutorials
+\brief Examples showing the TUnfold usage.
diff --git a/tutorials/unfold/testUnfold1.C b/tutorials/unfold/testUnfold1.C
new file mode 100644
index 00000000000..d32d2af1055
--- /dev/null
+++ b/tutorials/unfold/testUnfold1.C
@@ -0,0 +1,517 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfold and related.
+///
+/// 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
+///
+/// Systematic errors from the MC shape variation are included
+/// and propagated to the result
+///
+/// 3. fit the unfolded distribution, including the correlation matrix
+///
+/// 4. save six plots to a file `testUnfold1.ps`
+/// ~~~
+/// 1 2 3
+/// 4 5 6
+/// ~~~
+/// 1. 2d-plot of the matrix describing the migrations
+/// 2. generator-level distributions
+/// - blue: unfolded data, total errors
+/// - green: unfolded data, statistical errors
+/// - red: generated data
+/// - black: fit to green data points
+/// 3. detector level distributions
+/// - blue: unfolded data, folded back through the matrix
+/// - black: Monte Carlo (with wrong peal position)
+/// - blue: data
+/// 4. global correlation coefficients
+/// 5. \f$ \chi^2 \f$ as a function of \f$ log(\tau) \f$
+/// the star indicates the final choice of \f$ \tau \f$
+/// 6. the L curve
+///
+/// \macro_output
+/// \macro_image
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0, updated for using the classes TUnfoldDensity, TUnfoldBinning
+/// - Version 16.1, parallel to changes in TUnfold
+/// - Version 16.0, parallel to changes in TUnfold
+/// - Version 15, with automated L-curve scan
+/// - Version 14, with changes in TUnfoldSys.cxx
+/// - Version 13, include test of systematic errors
+/// - Version 12, catch error when defining the input
+/// - 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 and TUnfold.h
+/// - Version 8, with bug-fix in TUnfold.cxx and TUnfold.h
+/// - Version 7, with bug-fix in TUnfold.cxx and TUnfold.h
+/// - Version 6a, fix problem with dynamic array allocation under windows
+/// - Version 6, bug-fixes in TUnfold.C
+/// - Version 5, replace main() by testUnfold1()
+/// - 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+
+#include <TError.h>
+#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"
+
+// #define VERBOSE_LCURVE_SCAN
+
+using namespace std;
+
+TRandom *rnd=0;
+
+TH2 *gHistInvEMatrix;
+
+TVirtualFitter *gFitter=0;
+
+void chisquare_corr(Int_t &npar, Double_t * /*gin */, Double_t &f, Double_t *u, Int_t /*flag */) {
+ // Minimization function for H1s using a Chisquare method
+ // only one-dimensional histograms are supported
+ // Correlated errors are taken from an external inverse covariance matrix
+ // stored in a 2-dimensional histogram
+ Double_t x;
+ TH1 *hfit = (TH1*)gFitter->GetObjectFit();
+ TF1 *f1 = (TF1*)gFitter->GetUserFunc();
+
+ f1->InitArgs(&x,u);
+ npar = f1->GetNpar();
+ f = 0;
+
+ Int_t npfit = 0;
+ Int_t nPoints=hfit->GetNbinsX();
+ Double_t *df=new Double_t[nPoints];
+ for (Int_t i=0;i<nPoints;i++) {
+ x = hfit->GetBinCenter(i+1);
+ TF1::RejectPoint(kFALSE);
+ df[i] = f1->EvalPar(&x,u)-hfit->GetBinContent(i+1);
+ if (TF1::RejectedPoint()) df[i]=0.0;
+ else npfit++;
+ }
+ for (Int_t i=0;i<nPoints;i++) {
+ for (Int_t j=0;j<nPoints;j++) {
+ f += df[i]*df[j]*gHistInvEMatrix->GetBinContent(i+1,j+1);
+ }
+ }
+ delete[] df;
+ f1->SetNumberFitPoints(npfit);
+}
+
+Double_t bw_func(Double_t *x,Double_t *par) {
+ Double_t dm=x[0]-par[1];
+ return par[0]/(dm*dm+par[2]*par[2]);
+}
+
+
+// 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 testUnfold1()
+{
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ // show fit result
+ gStyle->SetOptFit(1111);
+
+ // 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 alternative MC
+ // this will be used to derive a systematic error due to MC
+ // parameter uncertainties
+ TH2D *histMdetGenSysMC=new TH2D("MdetgenSysMC",";mass(det);mass(gen)",
+ nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
+ neventMC=rnd->Poisson(luminosityMC*crossSection);
+ for(Int_t i=0;i<neventMC;i++) {
+ Double_t mGen=GenerateEvent
+ (0.5, // relative fraction of background
+ 3.6, // peak position in MC with systematic shift
+ 0.15); // peak width in MC
+ Double_t mDet=DetectorEvent(TMath::Abs(mGen));
+ histMdetGenSysMC->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 in data
+ 0.15); // peak width in data
+ 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);
+ }
+
+ //=========================================================================
+ // divide by bin width to get density distributions
+ TH1D *histDensityGenData=new TH1D("DensityGenData",";mass(gen)",
+ nGen,xminGen,xmaxGen);
+ TH1D *histDensityGenMC=new TH1D("DensityGenMC",";mass(gen)",
+ nGen,xminGen,xmaxGen);
+ for(Int_t i=1;i<=nGen;i++) {
+ histDensityGenData->SetBinContent(i,histMgenData->GetBinContent(i)/
+ histMgenData->GetBinWidth(i));
+ histDensityGenMC->SetBinContent(i,histMgenMC->GetBinContent(i)/
+ histMgenMC->GetBinWidth(i));
+ }
+
+ //=========================================================================
+ // set up the unfolding
+ // define migration matrix
+ TUnfoldDensity unfold(histMdetGenMC,TUnfold::kHistMapOutputVert);
+
+ // define input and bias scheme
+ // do not use the bias, because MC peak may be at the wrong place
+ // watch out for error codes returned by the SetInput method
+ // errors larger or equal 10000 are fatal:
+ // the data points specified as input are not sufficient to constrain the
+ // unfolding process
+ if(unfold.SetInput(histMdetData)>=10000) {
+ std::cout<<"Unfolding result may be wrong\n";
+ }
+
+ //========================================================================
+ // the unfolding is done here
+ //
+ // scan L curve and find best point
+ Int_t nScan=30;
+ // use automatic L-curve scan: start with taumin=taumax=0.0
+ Double_t tauMin=0.0;
+ Double_t tauMax=0.0;
+ Int_t iBest;
+ TSpline *logTauX,*logTauY;
+ TGraph *lCurve;
+
+ // if required, report Info messages (for debugging the L-curve scan)
+#ifdef VERBOSE_LCURVE_SCAN
+ Int_t oldinfo=gErrorIgnoreLevel;
+ gErrorIgnoreLevel=kInfo;
+#endif
+ // 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);
+
+ // if required, switch to previous log-level
+#ifdef VERBOSE_LCURVE_SCAN
+ gErrorIgnoreLevel=oldinfo;
+#endif
+
+ //==========================================================================
+ // define a correlated systematic error
+ // for example, assume there is a 10% correlated error for all reconstructed
+ // masses larger than 7
+ Double_t SYS_ERROR1_MSTART=6;
+ Double_t SYS_ERROR1_SIZE=0.1;
+ TH2D *histMdetGenSys1=new TH2D("Mdetgensys1",";mass(det);mass(gen)",
+ nDet,xminDet,xmaxDet,nGen,xminGen,xmaxGen);
+ for(Int_t i=0;i<=nDet+1;i++) {
+ if(histMdetData->GetBinCenter(i)>=SYS_ERROR1_MSTART) {
+ for(Int_t j=0;j<=nGen+1;j++) {
+ histMdetGenSys1->SetBinContent(i,j,SYS_ERROR1_SIZE);
+ }
+ }
+ }
+ unfold.AddSysError(histMdetGenSysMC,"SYSERROR_MC",TUnfold::kHistMapOutputVert,
+ TUnfoldSys::kSysErrModeMatrix);
+ unfold.AddSysError(histMdetGenSys1,"SYSERROR1",TUnfold::kHistMapOutputVert,
+ TUnfoldSys::kSysErrModeRelative);
+
+ //==========================================================================
+ // print some results
+ //
+ std::cout<<"tau="<<unfold.GetTau()<<"\n";
+ std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
+ <<" / "<<unfold.GetNdf()<<"\n";
+ std::cout<<"chi**2(sys)="<<unfold.GetChi2Sys()<<"\n";
+
+
+ //==========================================================================
+ // create graphs with one point to visualize the 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);
+
+ //==========================================================================
+ // retrieve results into histograms
+
+ // get unfolded distribution
+ TH1 *histMunfold=unfold.GetOutput("Unfolded");
+
+ // get unfolding result, folded back
+ TH1 *histMdetFold=unfold.GetFoldedOutput("FoldedBack");
+
+ // get error matrix (input distribution [stat] errors only)
+ // TH2D *histEmatData=unfold.GetEmatrix("EmatData");
+
+ // get total error matrix:
+ // migration matrix uncorrelated and correlated systematic errors
+ // added in quadrature to the data statistical errors
+ TH2 *histEmatTotal=unfold.GetEmatrixTotal("EmatTotal");
+
+ // create data histogram with the total errors
+ TH1D *histTotalError=
+ new TH1D("TotalError",";mass(gen)",nGen,xminGen,xmaxGen);
+ for(Int_t bin=1;bin<=nGen;bin++) {
+ histTotalError->SetBinContent(bin,histMunfold->GetBinContent(bin));
+ histTotalError->SetBinError
+ (bin,TMath::Sqrt(histEmatTotal->GetBinContent(bin,bin)));
+ }
+
+ // get global correlation coefficients
+ // for this calculation one has to specify whether the
+ // underflow/overflow bins are included or not
+ // default: include all bins
+ // here: exclude underflow and overflow bins
+ TH1 *histRhoi=unfold.GetRhoItotal("rho_I",
+ 0, // use default title
+ 0, // all distributions
+ "*[UO]", // discard underflow and overflow bins on all axes
+ kTRUE, // use original binning
+ &gHistInvEMatrix // store inverse of error matrix
+ );
+
+ //======================================================================
+ // fit Breit-Wigner shape to unfolded data, using the full error matrix
+ // here we use a "user" chi**2 function to take into account
+ // the full covariance matrix
+
+ gFitter=TVirtualFitter::Fitter(histMunfold);
+ gFitter->SetFCN(chisquare_corr);
+
+ TF1 *bw=new TF1("bw",bw_func,xminGen,xmaxGen,3);
+ bw->SetParameter(0,1000.);
+ bw->SetParameter(1,3.8);
+ bw->SetParameter(2,0.2);
+
+ // for (wrong!) fitting without correlations, drop the option "U"
+ // here.
+ histMunfold->Fit(bw,"UE");
+
+ //=====================================================================
+ // plot some histograms
+ TCanvas output;
+ 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);
+ histTotalError->SetLineColor(kBlue);
+ histTotalError->Draw("E");
+ histMunfold->SetLineColor(kGreen);
+ histMunfold->Draw("SAME E1");
+ histDensityGenData->SetLineColor(kRed);
+ histDensityGenData->Draw("SAME");
+ histDensityGenMC->Draw("SAME HIST");
+
+ // show detector level distributions
+ // data (red)
+ // MC (black) [with completely wrong peak position and shape]
+ // unfolded data (blue)
+ output.cd(3);
+ histMdetFold->SetLineColor(kBlue);
+ histMdetFold->Draw();
+ histMdetMC->Draw("SAME HIST");
+
+ TH1 *histInput=unfold.GetInput("Minput",";mass(det)");
+
+ histInput->SetLineColor(kRed);
+ histInput->Draw("SAME");
+
+ // show correlation coefficients
+ output.cd(4);
+ histRhoi->Draw();
+
+ // show tau as a function of chi**2
+ output.cd(5);
+ logTauX->Draw();
+ bestLogTauLogChi2->SetMarkerColor(kRed);
+ bestLogTauLogChi2->Draw("*");
+
+ // show the L curve
+ output.cd(6);
+ lCurve->Draw("AL");
+ bestLcurve->SetMarkerColor(kRed);
+ bestLcurve->Draw("*");
+
+ output.SaveAs("testUnfold1.ps");
+
+ return 0;
+}
diff --git a/tutorials/unfold/testUnfold2.C b/tutorials/unfold/testUnfold2.C
new file mode 100644
index 00000000000..a82fbbb747d
--- /dev/null
+++ b/tutorials/unfold/testUnfold2.C
@@ -0,0 +1,366 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// 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
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0, updated for changed methods in TUnfold
+/// - 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 separate 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#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;
+
+
+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 corresponding 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;
+}
diff --git a/tutorials/unfold/testUnfold3.C b/tutorials/unfold/testUnfold3.C
new file mode 100644
index 00000000000..bb83f95410e
--- /dev/null
+++ b/tutorials/unfold/testUnfold3.C
@@ -0,0 +1,631 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// 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
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0, change to use the class TUnfoldDensity
+/// - Version 16.1, parallel to changes in TUnfold
+/// - Version 16.0, parallel to changes in TUnfold
+/// - Version 15, simple example including background subtraction
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#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;
+
+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: retrieve 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)");
+
+ // retrieve error matrix of statistical errors
+ TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");
+
+ // retrieve 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);
+ }
+ }
+
+ // retrieve 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";
+ }
+}
diff --git a/tutorials/unfold/testUnfold4.C b/tutorials/unfold/testUnfold4.C
new file mode 100644
index 00000000000..f8a9a2a418b
--- /dev/null
+++ b/tutorials/unfold/testUnfold4.C
@@ -0,0 +1,222 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the class TUnfoldSys.
+///
+/// Simple toy tests of the TUnfold package
+///
+/// Pseudo data (5000 events) are unfolded into three components
+/// The unfolding is performed once without and once with area constraint
+///
+/// Ideally, the pulls may show that the result is biased if no constraint
+/// is applied. This is expected because the true data errors are not known,
+/// and instead the sqrt(data) errors are used.
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 16.1, parallel to changes in TUnfold
+/// - Version 16.0, parallel to changes in TUnfold
+/// - Version 15, use L-curve scan to scan the average correlation
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#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;
+
+TRandom *rnd=0;
+
+Int_t GenerateGenEvent(Int_t nmax,const Double_t *probability) {
+ // choose an integer random number in the range [0,nmax]
+ // (the generator level bin)
+ // depending on the probabilities
+ // probability[0],probability[1],...probability[nmax-1]
+ Double_t f=rnd->Rndm();
+ Int_t r=0;
+ while((r<nmax)&&(f>=probability[r])) {
+ f -= probability[r];
+ r++;
+ }
+ return r;
+}
+
+Double_t GenerateRecEvent(const Double_t *shapeParm) {
+ // return a coordinate (the reconstructed variable)
+ // depending on shapeParm[]
+ // shapeParm[0]: fraction of events with Gaussian distribution
+ // shapeParm[1]: mean of Gaussian
+ // shapeParm[2]: width of Gaussian
+ // (1-shapeParm[0]): fraction of events with flat distribution
+ // shapeParm[3]: minimum of flat component
+ // shapeParm[4]: maximum of flat component
+ Double_t f=rnd->Rndm();
+ Double_t r;
+ if(f<shapeParm[0]) {
+ r=rnd->Gaus(shapeParm[1],shapeParm[2]);
+ } else {
+ r=rnd->Rndm()*(shapeParm[4]-shapeParm[3])+shapeParm[3];
+ }
+ return r;
+}
+
+void testUnfold4()
+{
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ // random generator
+ rnd=new TRandom3();
+
+ // data and MC number of events
+ Double_t const nData0= 500.0;
+ Double_t const nMC0 = 50000.0;
+
+ // Binning
+ // reconstructed variable (0-10)
+ Int_t const nDet=15;
+ Double_t const xminDet=0.0;
+ Double_t const xmaxDet=15.0;
+
+ // signal binning (three shapes: 0,1,2)
+ Int_t const nGen=3;
+ Double_t const xminGen=-0.5;
+ Double_t const xmaxGen= 2.5;
+
+ // parameters
+ // fraction of events per signal shape
+ static const Double_t genFrac[]={0.3,0.6,0.1};
+
+ // signal shapes
+ static const Double_t genShape[][5]=
+ {{1.0,2.0,1.5,0.,15.},
+ {1.0,7.0,2.5,0.,15.},
+ {0.0,0.0,0.0,0.,15.}};
+
+ // define DATA histograms
+ // observed data distribution
+ TH1D *histDetDATA=new TH1D("Yrec",";DATA(Yrec)",nDet,xminDet,xmaxDet);
+
+ // define MC histograms
+ // matrix of migrations
+ TH2D *histGenDetMC=new TH2D("Yrec%Xgen","MC(Xgen,Yrec)",
+ nGen,xminGen,xmaxGen,nDet,xminDet,xmaxDet);
+
+ TH1D *histUnfold=new TH1D("Xgen",";DATA(Xgen)",nGen,xminGen,xmaxGen);
+
+ TH1D **histPullNC=new TH1D* [nGen];
+ TH1D **histPullArea=new TH1D* [nGen];
+ for(int i=0;i<nGen;i++) {
+ histPullNC[i]=new TH1D(TString::Format("PullNC%d",i),"pull",15,-3.,3.);
+ histPullArea[i]=new TH1D(TString::Format("PullArea%d",i),"pull",15,-3.,3.);
+ }
+
+ // this method is new in version 16 of TUnfold
+ cout<<"TUnfold version is "<<TUnfold::GetTUnfoldVersion()<<"\n";
+
+ for(int itoy=0;itoy<1000;itoy++) {
+ if(!(itoy %10)) cout<<"toy iteration: "<<itoy<<"\n";
+ histDetDATA->Reset();
+ histGenDetMC->Reset();
+
+ Int_t nData=rnd->Poisson(nData0);
+ for(Int_t i=0;i<nData;i++) {
+ Int_t iGen=GenerateGenEvent(nGen,genFrac);
+ Double_t yObs=GenerateRecEvent(genShape[iGen]);
+ histDetDATA->Fill(yObs);
+ }
+
+ Int_t nMC=rnd->Poisson(nMC0);
+ for(Int_t i=0;i<nMC;i++) {
+ Int_t iGen=GenerateGenEvent(nGen,genFrac);
+ Double_t yObs=GenerateRecEvent(genShape[iGen]);
+ histGenDetMC->Fill(iGen,yObs);
+ }
+ /* for(Int_t ix=0;ix<=histGenDetMC->GetNbinsX()+1;ix++) {
+ for(Int_t iy=0;iy<=histGenDetMC->GetNbinsY()+1;iy++) {
+ cout<<ix<<iy<<" : "<<histGenDetMC->GetBinContent(ix,iy)<<"\n";
+ }
+ } */
+ //========================
+ // unfolding
+
+ TUnfoldSys unfold(histGenDetMC,TUnfold::kHistMapOutputHoriz,
+ TUnfold::kRegModeSize,TUnfold::kEConstraintNone);
+ // define the input vector (the measured data distribution)
+ unfold.SetInput(histDetDATA,0.0,1.0);
+
+ // run the unfolding
+ unfold.ScanLcurve(50,0.,0.,0,0,0);
+
+ // fill pull distributions without constraint
+ unfold.GetOutput(histUnfold);
+
+ for(int i=0;i<nGen;i++) {
+ histPullNC[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
+ histUnfold->GetBinError(i+1));
+ }
+
+ // repeat unfolding on the same data, now with Area constraint
+ unfold.SetConstraint(TUnfold::kEConstraintArea);
+
+ // run the unfolding
+ unfold.ScanLcurve(50,0.,0.,0,0,0);
+
+ // fill pull distributions with constraint
+ unfold.GetOutput(histUnfold);
+
+ for(int i=0;i<nGen;i++) {
+ histPullArea[i]->Fill((histUnfold->GetBinContent(i+1)-genFrac[i]*nData0)/
+ histUnfold->GetBinError(i+1));
+ }
+
+ }
+ TCanvas output;
+ output.Divide(3,2);
+
+ gStyle->SetOptFit(1111);
+
+ for(int i=0;i<nGen;i++) {
+ output.cd(i+1);
+ histPullNC[i]->Fit("gaus");
+ histPullNC[i]->Draw();
+ }
+ for(int i=0;i<nGen;i++) {
+ output.cd(i+4);
+ histPullArea[i]->Fit("gaus");
+ histPullArea[i]->Draw();
+ }
+ output.SaveAs("testUnfold4.ps");
+}
diff --git a/tutorials/unfold/testUnfold5a.C b/tutorials/unfold/testUnfold5a.C
new file mode 100644
index 00000000000..ce57172e0bc
--- /dev/null
+++ b/tutorials/unfold/testUnfold5a.C
@@ -0,0 +1,338 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning
+///
+/// A toy test of the TUnfold package
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// - testUnfold5a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - testUnfold5b.C create a root file with the TUnfoldBinning objects
+/// - write file testUnfold5_binning.root
+///
+/// - testUnfold5c.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold5_binning.root
+/// testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - write testUnfold5_histograms.root
+///
+/// - testUnfold5d.C run the unfolding
+/// - read testUnfold5_histograms.root
+/// - write testUnfold5_result.root
+/// testUnfold5_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0 example for multi-dimensional unfolding
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <map>
+#include <cmath>
+#include <TMath.h>
+#include <TRandom3.h>
+#include <TFile.h>
+#include <TTree.h>
+
+using namespace std;
+
+TRandom *g_rnd=0;
+
+class ToyEvent {
+public:
+ void GenerateDataEvent(TRandom *rnd);
+ void GenerateSignalEvent(TRandom *rnd);
+ void GenerateBgrEvent(TRandom *rnd);
+ // reconstructed quantities
+ inline Double_t GetPtRec(void) const { return fPtRec; }
+ inline Double_t GetEtaRec(void) const { return fEtaRec; }
+ inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
+ inline Bool_t IsTriggered(void) const { return fIsTriggered; }
+
+ // generator level quantities
+ inline Double_t GetPtGen(void) const {
+ if(IsSignal()) return fPtGen;
+ else return -1.0;
+ }
+ inline Double_t GetEtaGen(void) const {
+ if(IsSignal()) return fEtaGen;
+ else return 999.0;
+ }
+ inline Bool_t IsSignal(void) const { return fIsSignal; }
+protected:
+
+ void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
+ void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
+ void GenerateReco(TRandom *rnd);
+
+ // reconstructed quantities
+ Double_t fPtRec;
+ Double_t fEtaRec;
+ Double_t fDiscriminator;
+ Bool_t fIsTriggered;
+ // generated quantities
+ Double_t fPtGen;
+ Double_t fEtaGen;
+ Bool_t fIsSignal;
+
+ static Double_t kDataSignalFraction;
+
+};
+
+void testUnfold5a()
+{
+ // random generator
+ g_rnd=new TRandom3();
+
+ // data and MC number of events
+ const Int_t neventData = 20000;
+ Double_t const neventSignalMC =2000000;
+ Double_t const neventBgrMC =2000000;
+
+ Float_t etaRec,ptRec,discr,etaGen,ptGen;
+ Int_t istriggered,issignal;
+
+ //==================================================================
+ // Step 1: generate data TTree
+
+ TFile *dataFile=new TFile("testUnfold5_data.root","recreate");
+ TTree *dataTree=new TTree("data","event");
+
+ dataTree->Branch("etarec",&etaRec,"etarec/F");
+ dataTree->Branch("ptrec",&ptRec,"ptrec/F");
+ dataTree->Branch("discr",&discr,"discr/F");
+
+ // for real data, only the triggered events are available
+ dataTree->Branch("istriggered",&istriggered,"istriggered/I");
+ // data truth parameters
+ dataTree->Branch("etagen",&etaGen,"etagen/F");
+ dataTree->Branch("ptgen",&ptGen,"ptgen/F");
+ dataTree->Branch("issignal",&issignal,"issignal/I");
+
+ cout<<"fill data tree\n";
+
+ Int_t nEvent=0,nTriggered=0;
+ while(nTriggered<neventData) {
+ ToyEvent event;
+ event.GenerateDataEvent(g_rnd);
+
+ etaRec=event.GetEtaRec();
+ ptRec=event.GetPtRec();
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+ etaGen=event.GetEtaGen();
+ ptGen=event.GetPtGen();
+ issignal=event.IsSignal() ? 1 : 0;
+
+ dataTree->Fill();
+
+ if(!(nEvent%100000)) cout<<" data event "<<nEvent<<"\n";
+
+ if(istriggered) nTriggered++;
+ nEvent++;
+
+ }
+
+ dataTree->Write();
+ delete dataTree;
+ delete dataFile;
+
+ //==================================================================
+ // Step 2: generate signal TTree
+
+ TFile *signalFile=new TFile("testUnfold5_signal.root","recreate");
+ TTree *signalTree=new TTree("signal","event");
+
+ signalTree->Branch("etarec",&etaRec,"etarec/F");
+ signalTree->Branch("ptrec",&ptRec,"ptrec/F");
+ signalTree->Branch("discr",&discr,"discr/F");
+ signalTree->Branch("istriggered",&istriggered,"istriggered/I");
+ signalTree->Branch("etagen",&etaGen,"etagen/F");
+ signalTree->Branch("ptgen",&ptGen,"ptgen/F");
+
+ cout<<"fill signal tree\n";
+
+ for(int ievent=0;ievent<neventSignalMC;ievent++) {
+ ToyEvent event;
+ event.GenerateSignalEvent(g_rnd);
+
+ etaRec=event.GetEtaRec();
+ ptRec=event.GetPtRec();
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+ etaGen=event.GetEtaGen();
+ ptGen=event.GetPtGen();
+
+ if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
+
+ signalTree->Fill();
+ }
+
+ signalTree->Write();
+ delete signalTree;
+ delete signalFile;
+
+ // ==============================================================
+ // Step 3: generate background MC TTree
+
+ TFile *bgrFile=new TFile("testUnfold5_background.root","recreate");
+ TTree *bgrTree=new TTree("background","event");
+
+ bgrTree->Branch("etarec",&etaRec,"etarec/F");
+ bgrTree->Branch("ptrec",&ptRec,"ptrec/F");
+ bgrTree->Branch("discr",&discr,"discr/F");
+ bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
+
+ cout<<"fill background tree\n";
+
+ for(int ievent=0;ievent<neventBgrMC;ievent++) {
+ ToyEvent event;
+ event.GenerateBgrEvent(g_rnd);
+ etaRec=event.GetEtaRec();
+ ptRec=event.GetPtRec();
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+
+ if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
+
+ bgrTree->Fill();
+ }
+
+ bgrTree->Write();
+ delete bgrTree;
+ delete bgrFile;
+}
+
+Double_t ToyEvent::kDataSignalFraction=0.8;
+
+void ToyEvent::GenerateDataEvent(TRandom *rnd) {
+ fIsSignal=rnd->Uniform()<kDataSignalFraction;
+ if(IsSignal()) {
+ GenerateSignalKinematics(rnd,kTRUE);
+ } else {
+ GenerateBgrKinematics(rnd,kTRUE);
+ }
+ GenerateReco(rnd);
+}
+
+void ToyEvent::GenerateSignalEvent(TRandom *rnd) {
+ fIsSignal=1;
+ GenerateSignalKinematics(rnd,kFALSE);
+ GenerateReco(rnd);
+}
+
+void ToyEvent::GenerateBgrEvent(TRandom *rnd) {
+ fIsSignal=0;
+ GenerateBgrKinematics(rnd,kFALSE);
+ GenerateReco(rnd);
+}
+
+void ToyEvent::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
+ Double_t e_T0=0.5;
+ Double_t e_T0_eta=0.0;
+ Double_t e_n=2.0;
+ Double_t e_n_eta=0.0;
+ Double_t eta_p2=0.0;
+ Double_t etaMax=4.0;
+ if(isData) {
+ e_T0=0.6;
+ e_n=2.5;
+ e_T0_eta=0.05;
+ e_n_eta=-0.05;
+ eta_p2=1.5;
+ }
+ if(eta_p2>0.0) {
+ fEtaGen=TMath::Power(rnd->Uniform(),eta_p2)*etaMax;
+ if(rnd->Uniform()>=0.5) fEtaGen= -fEtaGen;
+ } else {
+ fEtaGen=rnd->Uniform(-etaMax,etaMax);
+ }
+ Double_t n=e_n + e_n_eta*fEtaGen;
+ Double_t T0=e_T0 + e_T0_eta*fEtaGen;
+ fPtGen=(TMath::Power(rnd->Uniform(),1./(1.-n))-1.)*T0;
+ /* static int print=100;
+ if(print) {
+ cout<<fEtaGen
+ <<" "<<fPtGen
+ <<"\n";
+ print--;
+ } */
+}
+
+void ToyEvent::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
+ fPtGen=0.0;
+ fEtaGen=0.0;
+ fPtRec=rnd->Exp(isData ? 2.5 : 2.5);
+ fEtaRec=rnd->Uniform(-3.,3.);
+}
+
+void ToyEvent::GenerateReco(TRandom *rnd) {
+ if(fIsSignal) {
+ Double_t expEta=TMath::Exp(fEtaGen);
+ Double_t eGen=fPtGen*(expEta+1./expEta);
+ Double_t sigmaE=0.1*TMath::Sqrt(eGen)+(0.01+0.002*TMath::Abs(fEtaGen))
+ *eGen;
+ Double_t eRec;
+ do {
+ eRec=rnd->Gaus(eGen,sigmaE);
+ } while(eRec<=0.0);
+ Double_t sigmaEta=0.1+0.02*fEtaGen;
+ fEtaRec=rnd->Gaus(fEtaGen,sigmaEta);
+ fPtRec=eRec/(expEta+1./expEta);
+ do {
+ Double_t tauDiscr=0.08-0.04/(1.+fPtRec/10.0);
+ Double_t sigmaDiscr=0.01;
+ fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
+ } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
+ /* static int print=100;
+ if(print) {
+ cout<<fEtaGen<<" "<<fPtGen
+ <<" -> "<<fEtaRec<<" "<<fPtRec
+ <<"\n";
+ print--;
+ } */
+ } else {
+ do {
+ Double_t tauDiscr=0.15-0.05/(1.+fPtRec/5.0)+0.1*fEtaRec;
+ Double_t sigmaDiscr=0.02+0.01*fEtaRec;
+ fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
+ } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
+ }
+ fIsTriggered=(rnd->Uniform()<1./(TMath::Exp(-fPtRec+3.5)+1.));
+}
diff --git a/tutorials/unfold/testUnfold5b.C b/tutorials/unfold/testUnfold5b.C
new file mode 100644
index 00000000000..b96fc0791b0
--- /dev/null
+++ b/tutorials/unfold/testUnfold5b.C
@@ -0,0 +1,182 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
+///
+/// A toy test of the TUnfold package
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// - testUnfold5a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - testUnfold5b.C create a root file with the TUnfoldBinning objects
+/// - write file testUnfold5_binning.root
+///
+/// - testUnfold5c.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold5_binning.root
+/// testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - write testUnfold5_histograms.root
+///
+/// - testUnfold5d.C run the unfolding
+/// - read testUnfold5_histograms.root
+/// - write testUnfold5_result.root
+/// testUnfold5_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, updated for writing out XML code
+/// - Version 17.4, updated for writing out XML code
+/// - Version 17.3, updated for writing out XML code
+/// - Version 17.2, updated for writing out XML code
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0 example for multi-dimensional unfolding
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <fstream>
+#include <TFile.h>
+#include "TUnfoldBinningXML.h"
+#include <TF2.h>
+
+using namespace std;
+
+void testUnfold5b()
+{
+
+ // write binning schemes to root file
+ TFile *binningSchemes=new TFile("testUnfold5_binning.root","recreate");
+
+ // reconstructed pt, eta, discriminator
+#define NBIN_PT_FINE 8
+#define NBIN_ETA_FINE 10
+#define NBIN_DISCR 4
+
+ // generated pt, eta
+#define NBIN_PT_COARSE 3
+#define NBIN_ETA_COARSE 3
+
+ // pt binning
+ Double_t ptBinsFine[NBIN_PT_FINE+1]=
+ {3.5,4.0,4.5,5.0,6.0,7.0,8.0,10.0,13.0};
+ Double_t ptBinsCoarse[NBIN_PT_COARSE+1]=
+ { 4.0, 5.0, 7.0, 10.0};
+ // eta binning
+ Double_t etaBinsFine[NBIN_ETA_FINE+1]=
+ {-3.,-2.5,-2.0,-1.,-0.5,0.0,0.5,1.,2.,2.5,3.};
+ Double_t etaBinsCoarse[NBIN_ETA_COARSE+1]=
+ { -2.0, -0.5, 0.5, 2. };
+
+ // discriminator bins
+ Double_t discrBins[NBIN_DISCR+1]={0.,0.15,0.5,0.85,1.0};
+
+ //=======================================================================
+ // detector level binning scheme
+
+ TUnfoldBinning *detectorBinning=new TUnfoldBinning("detector");
+ // highest discriminator bin has fine binning
+ TUnfoldBinning *detectorDistribution=
+ detectorBinning->AddBinning("detectordistribution");
+ detectorDistribution->AddAxis("pt",NBIN_PT_FINE,ptBinsFine,
+ false, // no underflow bin (not reconstructed)
+ true // overflow bin
+ );
+ detectorDistribution->AddAxis("eta",NBIN_ETA_FINE,etaBinsFine,
+ false, // no underflow bin (not reconstructed)
+ false // no overflow bin (not reconstructed)
+ );
+ detectorDistribution->AddAxis("discriminator",NBIN_DISCR,discrBins,
+ false, // no underflow bin (empty)
+ false // no overflow bin (empty)
+ );
+ /* TUnfoldBinning *detectorExtra=
+ detectorBinning->AddBinning("detectorextra",7,"one;zwei;three");
+ detectorBinning->PrintStream(cout); */
+
+ //=======================================================================
+ // generator level binning
+ TUnfoldBinning *generatorBinning=new TUnfoldBinning("generator");
+
+ // signal distribution is measured with coarse binning
+ // underflow and overflow bins are needed ot take care of
+ // what happens outside the phase-space
+ TUnfoldBinning *signalBinning = generatorBinning->AddBinning("signal");
+ signalBinning->AddAxis("ptgen",NBIN_PT_COARSE,ptBinsCoarse,
+ true, // underflow bin
+ true // overflow bin
+ );
+ signalBinning->AddAxis("etagen",NBIN_ETA_COARSE,etaBinsCoarse,
+ true, // underflow bin
+ true // overflow bin
+ );
+
+ // this is just an example how to set bin-dependent factors
+ // for the regularisation
+ TF2 *userFunc=new TF2("userfunc","1./x+0.2*y^2",ptBinsCoarse[0],
+ ptBinsCoarse[NBIN_PT_COARSE],
+ etaBinsCoarse[0],etaBinsCoarse[NBIN_ETA_COARSE]);
+ signalBinning->SetBinFactorFunction(1.0,userFunc);
+
+ // background distribution is unfolded with fine binning
+ // !!! in the reconstructed variable !!!
+ //
+ // This has the effect of "normalizing" the background in each
+ // pt,eta bin to the low discriminator values
+ // Only the shape of the discriminator in each (pt,eta) bin
+ // is taken from Monte Carlo
+ //
+ // This method has been applied e.g. in
+ // H1 Collaboration, "Prompt photons in Photoproduction"
+ // Eur.Phys.J. C66 (2010) 17
+ //
+ TUnfoldBinning *bgrBinning = generatorBinning->AddBinning("background");
+ bgrBinning->AddAxis("ptrec",NBIN_PT_FINE,ptBinsFine,
+ false, // no underflow bin (not reconstructed)
+ true // overflow bin
+ );
+ bgrBinning->AddAxis("etarec",NBIN_ETA_FINE,etaBinsFine,
+ false, // no underflow bin (not reconstructed)
+ false // no overflow bin (not reconstructed)
+ );
+ generatorBinning->PrintStream(cout);
+
+ detectorBinning->Write();
+ generatorBinning->Write();
+
+ ofstream xmlOut("testUnfold5binning.xml");
+ TUnfoldBinningXML::ExportXML(*detectorBinning,xmlOut,kTRUE,kFALSE);
+ TUnfoldBinningXML::ExportXML(*generatorBinning,xmlOut,kFALSE,kTRUE);
+ TUnfoldBinningXML::WriteDTD();
+ xmlOut.close();
+
+ delete binningSchemes;
+}
diff --git a/tutorials/unfold/testUnfold5c.C b/tutorials/unfold/testUnfold5c.C
new file mode 100644
index 00000000000..ee09cda40c2
--- /dev/null
+++ b/tutorials/unfold/testUnfold5c.C
@@ -0,0 +1,290 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
+///
+/// A toy test of the TUnfold package
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// - testUnfold5a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - testUnfold5b.C create a root file with the TUnfoldBinning objects
+/// - write file testUnfold5_binning.root
+///
+/// - testUnfold5c.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold5_binning.root
+/// testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// - write testUnfold5_histograms.root
+///
+/// - testUnfold5d.C run the unfolding
+/// - read testUnfold5_histograms.root
+/// - write testUnfold5_result.root
+/// testUnfold5_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, updated for reading binning from XML file
+/// - Version 17.4, updated for reading binning from XML file
+/// - Version 17.3, updated for reading binning from XML file
+/// - Version 17.2, updated for reading binning from XML file
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0 example for multi-dimensional unfolding
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+// uncomment this to read the binning schemes from the root file
+// by default the binning is read from the XML file
+// #define READ_BINNING_CINT
+
+
+#include <iostream>
+#include <map>
+#include <cmath>
+#include <TMath.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TH1.h>
+#ifndef READ_BINNING_CINT
+#include <TDOMParser.h>
+#include <TXMLDocument.h>
+#include "TUnfoldBinningXML.h"
+#else
+#include "TUnfoldBinning.h"
+#endif
+
+using namespace std;
+
+
+void testUnfold5c()
+{
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ //=======================================================
+ // Step 1: open file to save histograms and binning schemes
+
+ TFile *outputFile=new TFile("testUnfold5_histograms.root","recreate");
+
+ //=======================================================
+ // Step 2: read binning from XML
+ // and save them to output file
+
+#ifdef READ_BINNING_CINT
+ TFile *binningSchemes=new TFile("testUnfold5_binning.root");
+#endif
+
+ TUnfoldBinning *detectorBinning,*generatorBinning;
+
+ outputFile->cd();
+
+ // read binning schemes in XML format
+#ifndef READ_BINNING_CINT
+ TDOMParser parser;
+ Int_t error=parser.ParseFile("testUnfold5binning.xml");
+ if(error) cout<<"error="<<error<<" from TDOMParser\n";
+ TXMLDocument const *XMLdocument=parser.GetXMLDocument();
+ detectorBinning=
+ TUnfoldBinningXML::ImportXML(XMLdocument,"detector");
+ generatorBinning=
+ TUnfoldBinningXML::ImportXML(XMLdocument,"generator");
+#else
+ binningSchemes->GetObject("detector",detectorBinning);
+ binningSchemes->GetObject("generator",generatorBinning);
+
+ delete binningSchemes;
+#endif
+ detectorBinning->Write();
+ generatorBinning->Write();
+
+ if(detectorBinning) {
+ detectorBinning->PrintStream(cout);
+ } else {
+ cout<<"could not read 'detector' binning\n";
+ }
+ if(generatorBinning) {
+ generatorBinning->PrintStream(cout);
+ } else {
+ cout<<"could not read 'generator' binning\n";
+ }
+
+ // pointers to various nodes in the binning scheme
+ const TUnfoldBinning *detectordistribution=
+ detectorBinning->FindNode("detectordistribution");
+
+ const TUnfoldBinning *signalBinning=
+ generatorBinning->FindNode("signal");
+
+ const TUnfoldBinning *bgrBinning=
+ generatorBinning->FindNode("background");
+
+ // write binning schemes to output file
+
+ //=======================================================
+ // Step 3: book and fill data histograms
+
+ Float_t etaRec,ptRec,discr,etaGen,ptGen;
+ Int_t istriggered,issignal;
+
+ outputFile->cd();
+
+ TH1 *histDataReco=detectorBinning->CreateHistogram("histDataReco");
+ TH1 *histDataTruth=generatorBinning->CreateHistogram("histDataTruth");
+
+ TFile *dataFile=new TFile("testUnfold5_data.root");
+ TTree *dataTree=(TTree *) dataFile->Get("data");
+
+ if(!dataTree) {
+ cout<<"could not read 'data' tree\n";
+ }
+
+ dataTree->ResetBranchAddresses();
+ dataTree->SetBranchAddress("etarec",&etaRec);
+ dataTree->SetBranchAddress("ptrec",&ptRec);
+ dataTree->SetBranchAddress("discr",&discr);
+ // for real data, only the triggered events are available
+ dataTree->SetBranchAddress("istriggered",&istriggered);
+ // data truth parameters
+ dataTree->SetBranchAddress("etagen",&etaGen);
+ dataTree->SetBranchAddress("ptgen",&ptGen);
+ dataTree->SetBranchAddress("issignal",&issignal);
+ dataTree->SetBranchStatus("*",1);
+
+
+ cout<<"loop over data events\n";
+
+ for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
+ if(dataTree->GetEntry(ievent)<=0) break;
+ // fill histogram with reconstructed quantities
+ if(istriggered) {
+ Int_t binNumber=
+ detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
+ histDataReco->Fill(binNumber);
+ }
+ // fill histogram with data truth parameters
+ if(issignal) {
+ // signal has true eta and pt
+ Int_t binNumber=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
+ histDataTruth->Fill(binNumber);
+ } else {
+ // background only has reconstructed pt and eta
+ Int_t binNumber=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
+ histDataTruth->Fill(binNumber);
+ }
+ }
+
+ delete dataTree;
+ delete dataFile;
+
+ //=======================================================
+ // Step 4: book and fill histogram of migrations
+ // it receives events from both signal MC and background MC
+
+ outputFile->cd();
+
+ TH2 *histMCGenRec=TUnfoldBinning::CreateHistogramOfMigrations
+ (generatorBinning,detectorBinning,"histMCGenRec");
+
+ TFile *signalFile=new TFile("testUnfold5_signal.root");
+ TTree *signalTree=(TTree *) signalFile->Get("signal");
+
+ if(!signalTree) {
+ cout<<"could not read 'signal' tree\n";
+ }
+
+ signalTree->ResetBranchAddresses();
+ signalTree->SetBranchAddress("etarec",&etaRec);
+ signalTree->SetBranchAddress("ptrec",&ptRec);
+ signalTree->SetBranchAddress("discr",&discr);
+ signalTree->SetBranchAddress("istriggered",&istriggered);
+ signalTree->SetBranchAddress("etagen",&etaGen);
+ signalTree->SetBranchAddress("ptgen",&ptGen);
+ signalTree->SetBranchStatus("*",1);
+
+ cout<<"loop over MC signal events\n";
+
+ for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
+ if(signalTree->GetEntry(ievent)<=0) break;
+
+ // bin number on generator level for signal
+ Int_t genBin=signalBinning->GetGlobalBinNumber(ptGen,etaGen);
+
+ // bin number on reconstructed level
+ // bin number 0 corresponds to non-reconstructed events
+ Int_t recBin=0;
+ if(istriggered) {
+ recBin=detectordistribution->GetGlobalBinNumber(ptRec,etaRec,discr);
+ }
+ histMCGenRec->Fill(genBin,recBin);
+ }
+
+ delete signalTree;
+ delete signalFile;
+
+ TFile *bgrFile=new TFile("testUnfold5_background.root");
+ TTree *bgrTree=(TTree *) bgrFile->Get("background");
+
+ if(!bgrTree) {
+ cout<<"could not read 'background' tree\n";
+ }
+
+ bgrTree->ResetBranchAddresses();
+ bgrTree->SetBranchAddress("etarec",&etaRec);
+ bgrTree->SetBranchAddress("ptrec",&ptRec);
+ bgrTree->SetBranchAddress("discr",&discr);
+ bgrTree->SetBranchAddress("istriggered",&istriggered);
+ bgrTree->SetBranchStatus("*",1);
+
+ cout<<"loop over MC background events\n";
+
+ for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
+ if(bgrTree->GetEntry(ievent)<=0) break;
+
+ // here, for background only reconstructed quantities are known
+ // and only the reconstructed events are relevant
+ if(istriggered) {
+ // bin number on generator level for background
+ Int_t genBin=bgrBinning->GetGlobalBinNumber(ptRec,etaRec);
+ // bin number on reconstructed level
+ Int_t recBin=detectordistribution->GetGlobalBinNumber
+ (ptRec,etaRec,discr);
+ histMCGenRec->Fill(genBin,recBin);
+ }
+ }
+
+ delete bgrTree;
+ delete bgrFile;
+
+ outputFile->Write();
+ delete outputFile;
+
+}
diff --git a/tutorials/unfold/testUnfold5d.C b/tutorials/unfold/testUnfold5d.C
new file mode 100644
index 00000000000..fd310077708
--- /dev/null
+++ b/tutorials/unfold/testUnfold5d.C
@@ -0,0 +1,296 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
+///
+/// A toy test of the TUnfold package
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// testUnfold5a.C create root files with TTree objects for
+/// signal, background and data
+/// -> write files testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// testUnfold5b.C create a root file with the TUnfoldBinning objects
+/// -> write file testUnfold5_binning.root
+///
+/// testUnfold5c.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// -> read testUnfold5_binning.root
+/// testUnfold5_signal.root
+/// testUnfold5_background.root
+/// testUnfold5_data.root
+///
+/// -> write testUnfold5_histograms.root
+///
+/// testUnfold5d.C run the unfolding
+/// -> read testUnfold5_histograms.root
+/// -> write testUnfold5_result.root
+/// testUnfold5_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, in parallel to changes in TUnfold
+/// - Version 17.2, in parallel to changes in TUnfold
+/// - Version 17.1, in parallel to changes in TUnfold
+/// - Version 17.0 example for multi-dimensional unfolding
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <cmath>
+#include <map>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TGraph.h>
+#include <TFile.h>
+#include <TH1.h>
+#include "TUnfoldDensity.h"
+
+using namespace std;
+
+// #define PRINT_MATRIX_L
+
+#define TEST_INPUT_COVARIANCE
+
+void testUnfold5d()
+{
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ //==============================================
+ // step 1 : open output file
+ TFile *outputFile=new TFile("testUnfold5_results.root","recreate");
+
+ //==============================================
+ // step 2 : read binning schemes and input histograms
+ TFile *inputFile=new TFile("testUnfold5_histograms.root");
+
+ outputFile->cd();
+
+ TUnfoldBinning *detectorBinning,*generatorBinning;
+
+ inputFile->GetObject("detector",detectorBinning);
+ inputFile->GetObject("generator",generatorBinning);
+
+ if((!detectorBinning)||(!generatorBinning)) {
+ cout<<"problem to read binning schemes\n";
+ }
+
+ // save binning schemes to output file
+ detectorBinning->Write();
+ generatorBinning->Write();
+
+ // read histograms
+ TH1 *histDataReco,*histDataTruth;
+ TH2 *histMCGenRec;
+
+ inputFile->GetObject("histDataReco",histDataReco);
+ inputFile->GetObject("histDataTruth",histDataTruth);
+ inputFile->GetObject("histMCGenRec",histMCGenRec);
+
+#ifdef TEST_ZERO_UNCORR_ERROR
+ // special test (bug in version 17.2 and below)
+ // set all errors in hisMCGenRec to zero
+ // -> program will crash
+ for(int i=0;i<=histMCGenRec->GetNbinsX()+1;i++) {
+ for(int j=0;j<=histMCGenRec->GetNbinsY()+1;j++) {
+ histMCGenRec->SetBinError(i,j,0.0);
+ }
+ }
+#endif
+
+ histDataReco->Write();
+ histDataTruth->Write();
+ histMCGenRec->Write();
+
+ if((!histDataReco)||(!histDataTruth)||(!histMCGenRec)) {
+ cout<<"problem to read input histograms\n";
+ }
+
+ //========================
+ // Step 3: unfolding
+
+ // preserve the area
+ TUnfold::EConstraint constraintMode= TUnfold::kEConstraintArea;
+
+ // basic choice of regularisation scheme:
+ // curvature (second derivative)
+ TUnfold::ERegMode regMode = TUnfold::kRegModeCurvature;
+
+ // density flags
+ TUnfoldDensity::EDensityMode densityFlags=
+ TUnfoldDensity::kDensityModeBinWidth;
+
+ // detailed steering for regularisation
+ const char *REGULARISATION_DISTRIBUTION=0;
+ const char *REGULARISATION_AXISSTEERING="*[B]";
+
+ // set up matrix of migrations
+ TUnfoldDensity unfold(histMCGenRec,TUnfold::kHistMapOutputHoriz,
+ regMode,constraintMode,densityFlags,
+ generatorBinning,detectorBinning,
+ REGULARISATION_DISTRIBUTION,
+ REGULARISATION_AXISSTEERING);
+
+ // define the input vector (the measured data distribution)
+
+#ifdef TEST_INPUT_COVARIANCE
+ // special test to use input covariance matrix
+ TH2D *inputEmatrix=
+ detectorBinning->CreateErrorMatrixHistogram("input_covar",true);
+ for(int i=1;i<=inputEmatrix->GetNbinsX();i++) {
+ Double_t e=histDataReco->GetBinError(i);
+ inputEmatrix->SetBinContent(i,i,e*e);
+ // test: non-zero covariance where variance is zero
+ // if(e<=0.) inputEmatrix->SetBinContent(i,i+1,1.0);
+ }
+ unfold.SetInput(histDataReco,0.0,0.0,inputEmatrix);
+#else
+ unfold.SetInput(histDataReco /* ,0.0,1.0 */);
+#endif
+ // print matrix of regularisation conditions
+#ifdef PRINT_MATRIX_L
+ TH2 *histL= unfold.GetL("L");
+ for(Int_t j=1;j<=histL->GetNbinsY();j++) {
+ cout<<"L["<<unfold.GetLBinning()->GetBinName(j)<<"]";
+ for(Int_t i=1;i<=histL->GetNbinsX();i++) {
+ Double_t c=histL->GetBinContent(i,j);
+ if(c!=0.0) cout<<" ["<<i<<"]="<<c;
+ }
+ cout<<"\n";
+ }
+#endif
+ // run the unfolding
+ //
+ // here, tau is determined by scanning the global correlation coefficients
+
+ Int_t nScan=30;
+ TSpline *rhoLogTau=0;
+ TGraph *lCurve=0;
+
+ // for determining tau, scan the correlation coefficients
+ // correlation coefficients may be probed for all distributions
+ // or only for selected distributions
+ // underflow/overflow bins may be included/excluded
+ //
+ const char *SCAN_DISTRIBUTION="signal";
+ const char *SCAN_AXISSTEERING=0;
+
+ Int_t iBest=unfold.ScanTau(nScan,0.,0.,&rhoLogTau,
+ TUnfoldDensity::kEScanTauRhoMax,
+ SCAN_DISTRIBUTION,SCAN_AXISSTEERING,
+ &lCurve);
+
+ // create graphs with one point to visualize best choice of tau
+ Double_t t[1],rho[1],x[1],y[1];
+ rhoLogTau->GetKnot(iBest,t[0],rho[0]);
+ lCurve->GetPoint(iBest,x[0],y[0]);
+ TGraph *bestRhoLogTau=new TGraph(1,t,rho);
+ TGraph *bestLCurve=new TGraph(1,x,y);
+ Double_t *tAll=new Double_t[nScan],*rhoAll=new Double_t[nScan];
+ for(Int_t i=0;i<nScan;i++) {
+ rhoLogTau->GetKnot(i,tAll[i],rhoAll[i]);
+ }
+ TGraph *knots=new TGraph(nScan,tAll,rhoAll);
+
+ cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
+ <<" / "<<unfold.GetNdf()<<"\n";
+
+
+ //===========================
+ // Step 4: retrieve and plot unfolding results
+
+ // get unfolding output
+ TH1 *histDataUnfold=unfold.GetOutput("unfolded signal",0,0,0,kFALSE);
+ // get Monte Carlo reconstructed data
+ TH1 *histMCReco=histMCGenRec->ProjectionY("histMCReco",0,-1,"e");
+ TH1 *histMCTruth=histMCGenRec->ProjectionX("histMCTruth",0,-1,"e");
+ Double_t scaleFactor=histDataTruth->GetSumOfWeights()/
+ histMCTruth->GetSumOfWeights();
+ histMCReco->Scale(scaleFactor);
+ histMCTruth->Scale(scaleFactor);
+ // get matrix of probabilities
+ TH2 *histProbability=unfold.GetProbabilityMatrix("histProbability");
+ // get global correlation coefficients
+ /* TH1 *histGlobalCorr=*/ unfold.GetRhoItotal("histGlobalCorr",0,0,0,kFALSE);
+ TH1 *histGlobalCorrScan=unfold.GetRhoItotal
+ ("histGlobalCorrScan",0,SCAN_DISTRIBUTION,SCAN_AXISSTEERING,kFALSE);
+ /* TH2 *histCorrCoeff=*/ unfold.GetRhoIJtotal("histCorrCoeff",0,0,0,kFALSE);
+
+ TCanvas canvas;
+ canvas.Print("testUnfold5.ps[");
+
+ //========== page 1 ============
+ // unfolding control plots
+ // input, matrix, output
+ // tau-scan, global correlations, correlation coefficients
+ canvas.Clear();
+ canvas.Divide(3,2);
+
+ // (1) all bins, compare to original MC distribution
+ canvas.cd(1);
+ histDataReco->SetMinimum(0.0);
+ histDataReco->Draw("E");
+ histMCReco->SetLineColor(kBlue);
+ histMCReco->Draw("SAME HIST");
+ // (2) matrix of probabilities
+ canvas.cd(2);
+ histProbability->Draw("BOX");
+ // (3) unfolded data, data truth, MC truth
+ canvas.cd(3);
+ gPad->SetLogy();
+ histDataUnfold->Draw("E");
+ histDataTruth->SetLineColor(kBlue);
+ histDataTruth->Draw("SAME HIST");
+ histMCTruth->SetLineColor(kRed);
+ histMCTruth->Draw("SAME HIST");
+ // (4) scan of correlation vs tau
+ canvas.cd(4);
+ rhoLogTau->Draw();
+ knots->Draw("*");
+ bestRhoLogTau->SetMarkerColor(kRed);
+ bestRhoLogTau->Draw("*");
+ // (5) global correlation coefficients for the distributions
+ // used during the scan
+ canvas.cd(5);
+ //histCorrCoeff->Draw("BOX");
+ histGlobalCorrScan->Draw("HIST");
+ // (6) L-curve
+ canvas.cd(6);
+ lCurve->Draw("AL");
+ bestLCurve->SetMarkerColor(kRed);
+ bestLCurve->Draw("*");
+
+
+ canvas.Print("testUnfold5.ps");
+
+ canvas.Print("testUnfold5.ps]");
+
+}
diff --git a/tutorials/unfold/testUnfold6.C b/tutorials/unfold/testUnfold6.C
new file mode 100644
index 00000000000..455af206d1c
--- /dev/null
+++ b/tutorials/unfold/testUnfold6.C
@@ -0,0 +1,117 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the class TUnfoldBinning.
+///
+/// read a simple binning scheme and create/test bin maps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// #### History:
+/// - Version 17.5, in parallel to changes in TUnfold
+/// - Version 17.4, in parallel to changes in TUnfold
+/// - Version 17.3, test bin map functionality in TUnfoldBinning
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <fstream>
+#include <iomanip>
+#include <map>
+#include <TDOMParser.h>
+#include <TXMLDocument.h>
+#include "TUnfoldBinningXML.h"
+
+using namespace std;
+
+void PrintBinMap(TUnfoldBinning *binning,const char * where,
+ const Int_t *binMap);
+
+void testUnfold6()
+{
+ TDOMParser parser;
+ ofstream dtdFile("tunfoldbinning.dtd");
+ TUnfoldBinningXML::WriteDTD(dtdFile);
+ dtdFile.close();
+ Int_t error=parser.ParseFile("testUnfold6binning.xml");
+ if(error) cout<<"error="<<error<<" from TDOMParser\n";
+ TXMLDocument const *XMLdocument=parser.GetXMLDocument();
+ TUnfoldBinningXML *binning=
+ TUnfoldBinningXML::ImportXML(XMLdocument,"binning");
+ if(!binning) {
+ cout<<"error: can not read binning (document empty?)\n";
+ } else {
+ cout<<"Binning scheme:\n =================================\n";
+ binning->PrintStream(cout);
+ Int_t *binMap=binning->CreateEmptyBinMap();
+ PrintBinMap(binning,"CreateEmptyBinMap",binMap);
+
+ TUnfoldBinning const *branch1=binning->FindNode("branch1");
+ branch1->FillBinMap1D(binMap,"y[C]",2);
+ PrintBinMap(binning,"branch1->FillBinMap1D(...,\"y[C]\",...,2)",binMap);
+
+ delete binMap; binMap=binning->CreateEmptyBinMap();
+ TUnfoldBinning const *branch2=binning->FindNode("branch2");
+ branch2->FillBinMap1D(binMap,"x[C]",7);
+ PrintBinMap(binning,"branch2->FillBinMap1D(...,\"x[C]\",...,7)",binMap);
+
+ delete binMap; binMap=binning->CreateEmptyBinMap();
+ binning->FillBinMap1D(binMap,"y[C]",1);
+ PrintBinMap(binning,"binning->FillBinMap1D(...,\"y[C]\",...,1)",binMap);
+
+ binning->ExportXML("testUnfold6.out.xml");
+ }
+}
+
+void PrintBinMap(TUnfoldBinning *binning,const char * where,
+ const Int_t *binMap) {
+
+ cout<<"\n"<<where<<"\n=======================\n";
+ cout<<"global bin:";
+ for(int i=0;i<binning->GetEndBin()+1;i++) {
+ cout<<setw(3)<<i;
+ }
+ cout<<"\n";
+ cout<<"mapped to: ";
+ for(int i=0;i<binning->GetEndBin()+1;i++) {
+ cout<<setw(3)<<binMap[i];
+ }
+ cout<<"\n";
+ map<int,vector<int> > destBin;
+ for(int i=0;i<binning->GetEndBin()+1;i++) {
+ destBin[binMap[i]].push_back(i);
+ }
+ bool printed=false;
+ for(map<int,vector<int> >::const_iterator i=destBin.begin();i!=destBin.end();i++) {
+ if((*i).first>=0) {
+ if(!printed) {
+ cout<<"\ndest |contributing bins\n"
+ <<"=====+======================================\n";
+ printed=true;
+ }
+ for(size_t j=0;j<(*i).second.size();j++) {
+ cout<<setw(4)<<(*i).first<<" |";
+ cout<<setw(3)<<binning->GetBinName((*i).second[j])<<"\n";
+ }
+ cout<<"=====+======================================\n";
+ }
+ }
+}
diff --git a/tutorials/unfold/testUnfold6binning.xml b/tutorials/unfold/testUnfold6binning.xml
new file mode 100644
index 00000000000..fc7cfda6637
--- /dev/null
+++ b/tutorials/unfold/testUnfold6binning.xml
@@ -0,0 +1,27 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
+<TUnfoldBinning>
+<BinningNode name="binning" firstbin="1" factor="1">
+ <BinningNode name="branch1" firstbin="1" factor="1">
+ <Axis name="x" lowEdge="0">
+ <Bin repeat="4" width="1" center="0.5" />
+ <Axis name="y" lowEdge="0">
+ <Bin repeat="3" width="1" center="0.5" />
+ </Axis>
+ </Axis>
+ </BinningNode>
+ <BinningNode name="branch2" firstbin="13" factor="1">
+ <Axis name="x" lowEdge="0">
+ <Bin repeat="3" width="1" center="0.5" />
+ <Axis name="y" lowEdge="0">
+ <Bin repeat="2" width="1" center="0.5" />
+ </Axis>
+ </Axis>
+ </BinningNode>
+ <BinningNode name="branch3" firstbin="19" factor="1">
+ <Axis name="x" lowEdge="0">
+ <Bin repeat="2" width="1" center="0.5" />
+ </Axis>
+ </BinningNode>
+</BinningNode>
+</TUnfoldBinning>
diff --git a/tutorials/unfold/testUnfold7a.C b/tutorials/unfold/testUnfold7a.C
new file mode 100644
index 00000000000..3aafe5f06d2
--- /dev/null
+++ b/tutorials/unfold/testUnfold7a.C
@@ -0,0 +1,447 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
+///
+/// A toy test of the TUnfold package
+///
+/// This example is documented in conference proceedings:
+///
+/// arXiv:1611.01927
+/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// - testUnfold7a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - testUnfold7b.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold7binning.xml
+/// testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - write testUnfold7_histograms.root
+///
+/// - testUnfold7c.C run the unfolding
+/// - read testUnfold7_histograms.root
+/// - write testUnfold7_result.root
+/// testUnfold7_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <map>
+#include <cmath>
+#include <TMath.h>
+#include <TRandom3.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TLorentzVector.h>
+
+#define MASS1 0.511E-3
+
+using namespace std;
+
+TRandom *g_rnd=0;
+
+class ToyEvent7 {
+public:
+ void GenerateDataEvent(TRandom *rnd);
+ void GenerateSignalEvent(TRandom *rnd);
+ void GenerateBgrEvent(TRandom *rnd);
+ // reconstructed quantities
+ inline Double_t GetMRec(int i) const { return fMRec[i]; }
+ inline Double_t GetPtRec(int i) const { return fPtRec[i]; }
+ inline Double_t GetEtaRec(int i) const { return fEtaRec[i]; }
+ inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
+ inline Double_t GetPhiRec(int i) const { return fPhiRec[i]; }
+ inline Bool_t IsTriggered(void) const { return fIsTriggered; }
+
+ // generator level quantities
+ inline Double_t GetMGen(int i) const {
+ if(IsSignal()) return fMGen[i];
+ else return -1.0;
+ }
+ inline Double_t GetPtGen(int i) const {
+ if(IsSignal()) return fPtGen[i];
+ else return -1.0;
+ }
+ inline Double_t GetEtaGen(int i) const {
+ if(IsSignal()) return fEtaGen[i];
+ else return 999.0;
+ }
+ inline Double_t GetPhiGen(int i) const {
+ if(IsSignal()) return fPhiGen[i];
+ else return 999.0;
+ }
+ inline Bool_t IsSignal(void) const { return fIsSignal; }
+protected:
+
+ void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
+ void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
+ void GenerateReco(TRandom *rnd);
+
+ // reconstructed quantities
+ Double_t fMRec[3];
+ Double_t fPtRec[3];
+ Double_t fEtaRec[3];
+ Double_t fPhiRec[3];
+ Double_t fDiscriminator;
+ Bool_t fIsTriggered;
+ // generated quantities
+ Double_t fMGen[3];
+ Double_t fPtGen[3];
+ Double_t fEtaGen[3];
+ Double_t fPhiGen[3];
+ Bool_t fIsSignal;
+public:
+ static Double_t kDataSignalFraction;
+ static Double_t kMCSignalFraction;
+
+};
+
+void testUnfold7a()
+{
+ // random generator
+ g_rnd=new TRandom3(4711);
+
+ // data and MC number of events
+ Double_t muData0=5000.;
+ // luminosity error
+ Double_t muData=muData0*g_rnd->Gaus(1.0,0.03);
+ // stat error
+ Int_t neventData = g_rnd->Poisson( muData);
+
+ // generated number of MC events
+ Int_t neventSigmc = 250000;
+ Int_t neventBgrmc = 100000;
+
+ Float_t etaRec[3],ptRec[3],phiRec[3],mRec[3],discr;
+ Float_t etaGen[3],ptGen[3],phiGen[3],mGen[3];
+ Float_t weight;
+ Int_t istriggered,issignal;
+
+ //==================================================================
+ // Step 1: generate data TTree
+
+ TFile *dataFile=new TFile("testUnfold7_data.root","recreate");
+ TTree *dataTree=new TTree("data","event");
+
+ dataTree->Branch("etarec",etaRec,"etarec[3]/F");
+ dataTree->Branch("ptrec",ptRec,"ptrec[3]/F");
+ dataTree->Branch("phirec",phiRec,"phirec[3]/F");
+ dataTree->Branch("mrec",mRec,"mrec[3]/F");
+ dataTree->Branch("discr",&discr,"discr/F");
+
+ // for real data, only the triggered events are available
+ dataTree->Branch("istriggered",&istriggered,"istriggered/I");
+ // data truth parameters
+ dataTree->Branch("etagen",etaGen,"etagen[3]/F");
+ dataTree->Branch("ptgen",ptGen,"ptgen[3]/F");
+ dataTree->Branch("phigen",phiGen,"phigen[3]/F");
+ dataTree->Branch("mgen",mGen,"mgen[3]/F");
+ dataTree->Branch("issignal",&issignal,"issignal/I");
+
+ cout<<"fill data tree\n";
+
+ //Int_t nEvent=0,nTriggered=0;
+ for(int ievent=0;ievent<neventData;ievent++) {
+ ToyEvent7 event;
+ event.GenerateDataEvent(g_rnd);
+ for(int i=0;i<3;i++) {
+ etaRec[i]=event.GetEtaRec(i);
+ ptRec[i]=event.GetPtRec(i);
+ phiRec[i]=event.GetPhiRec(i);
+ mRec[i]=event.GetMRec(i);
+ etaGen[i]=event.GetEtaGen(i);
+ ptGen[i]=event.GetPtGen(i);
+ phiGen[i]=event.GetPhiGen(i);
+ mGen[i]=event.GetMGen(i);
+ }
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+ issignal=event.IsSignal() ? 1 : 0;
+
+ dataTree->Fill();
+
+ if(!(ievent%100000)) cout<<" data event "<<ievent<<"\n";
+
+ //if(istriggered) nTriggered++;
+ //nEvent++;
+
+ }
+
+ dataTree->Write();
+ delete dataTree;
+ delete dataFile;
+
+ //==================================================================
+ // Step 2: generate signal TTree
+
+ TFile *signalFile=new TFile("testUnfold7_signal.root","recreate");
+ TTree *signalTree=new TTree("signal","event");
+
+ signalTree->Branch("etarec",etaRec,"etarec[3]/F");
+ signalTree->Branch("ptrec",ptRec,"ptrec[3]/F");
+ signalTree->Branch("phirec",ptRec,"phirec[3]/F");
+ signalTree->Branch("mrec",mRec,"mrec[3]/F");
+ signalTree->Branch("discr",&discr,"discr/F");
+
+ // for real data, only the triggered events are available
+ signalTree->Branch("istriggered",&istriggered,"istriggered/I");
+ // data truth parameters
+ signalTree->Branch("etagen",etaGen,"etagen[3]/F");
+ signalTree->Branch("ptgen",ptGen,"ptgen[3]/F");
+ signalTree->Branch("phigen",phiGen,"phigen[3]/F");
+ signalTree->Branch("weight",&weight,"weight/F");
+ signalTree->Branch("mgen",mGen,"mgen[3]/F");
+
+ cout<<"fill signal tree\n";
+
+ weight=ToyEvent7::kMCSignalFraction*muData0/neventSigmc;
+
+ for(int ievent=0;ievent<neventSigmc;ievent++) {
+ ToyEvent7 event;
+ event.GenerateSignalEvent(g_rnd);
+
+ for(int i=0;i<3;i++) {
+ etaRec[i]=event.GetEtaRec(i);
+ ptRec[i]=event.GetPtRec(i);
+ phiRec[i]=event.GetPhiRec(i);
+ mRec[i]=event.GetMRec(i);
+ etaGen[i]=event.GetEtaGen(i);
+ ptGen[i]=event.GetPtGen(i);
+ phiGen[i]=event.GetPhiGen(i);
+ mGen[i]=event.GetMGen(i);
+ }
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+
+ if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
+
+ signalTree->Fill();
+ }
+
+ signalTree->Write();
+ delete signalTree;
+ delete signalFile;
+
+ // ==============================================================
+ // Step 3: generate background MC TTree
+
+ TFile *bgrFile=new TFile("testUnfold7_background.root","recreate");
+ TTree *bgrTree=new TTree("background","event");
+
+ bgrTree->Branch("etarec",&etaRec,"etarec[3]/F");
+ bgrTree->Branch("ptrec",&ptRec,"ptrec[3]/F");
+ bgrTree->Branch("phirec",&phiRec,"phirec[3]/F");
+ bgrTree->Branch("mrec",&mRec,"mrec[3]/F");
+ bgrTree->Branch("discr",&discr,"discr/F");
+ bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
+ signalTree->Branch("weight",&weight,"weight/F");
+
+ cout<<"fill background tree\n";
+
+ weight=(1.-ToyEvent7::kMCSignalFraction)*muData0/neventBgrmc;
+
+ for(int ievent=0;ievent<neventBgrmc;ievent++) {
+ ToyEvent7 event;
+ event.GenerateBgrEvent(g_rnd);
+ for(int i=0;i<3;i++) {
+ etaRec[i]=event.GetEtaRec(i);
+ ptRec[i]=event.GetPtRec(i);
+ phiRec[i]=event.GetPhiRec(i);
+ }
+ discr=event.GetDiscriminator();
+ istriggered=event.IsTriggered() ? 1 : 0;
+
+ if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
+
+ bgrTree->Fill();
+ }
+
+ bgrTree->Write();
+ delete bgrTree;
+ delete bgrFile;
+}
+
+Double_t ToyEvent7::kDataSignalFraction=0.75;
+Double_t ToyEvent7::kMCSignalFraction=0.75;
+
+void ToyEvent7::GenerateDataEvent(TRandom *rnd) {
+ fIsSignal=rnd->Uniform()<kDataSignalFraction;
+ if(IsSignal()) {
+ GenerateSignalKinematics(rnd,kTRUE);
+ } else {
+ GenerateBgrKinematics(rnd,kTRUE);
+ }
+ GenerateReco(rnd);
+}
+
+void ToyEvent7::GenerateSignalEvent(TRandom *rnd) {
+ fIsSignal=1;
+ GenerateSignalKinematics(rnd,kFALSE);
+ GenerateReco(rnd);
+}
+
+void ToyEvent7::GenerateBgrEvent(TRandom *rnd) {
+ fIsSignal=0;
+ GenerateBgrKinematics(rnd,kFALSE);
+ GenerateReco(rnd);
+}
+
+void ToyEvent7::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
+
+ // fake decay of Z0 to two fermions
+ double M0=91.1876;
+ double Gamma=2.4952;
+ // generated mass
+ do {
+ fMGen[2]=rnd->BreitWigner(M0,Gamma);
+ } while(fMGen[2]<=0.0);
+
+ double N_ETA=3.0;
+ double MU_PT=5.;
+ double SIGMA_PT=2.0;
+ double DECAY_A=0.2;
+ if(isData) {
+ //N_ETA=2.5;
+ MU_PT=6.;
+ SIGMA_PT=1.8;
+ //DECAY_A=0.5;
+ }
+ fEtaGen[2]=TMath::Power(rnd->Uniform(0,1.5),N_ETA);
+ if(rnd->Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
+ fPhiGen[2]=rnd->Uniform(-M_PI,M_PI);
+ do {
+ fPtGen[2]=rnd->Landau(MU_PT,SIGMA_PT);
+ } while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
+ //========================== decay
+ TLorentzVector sum;
+ sum.SetPtEtaPhiM(fPtGen[2],fEtaGen[2],fPhiGen[2],fMGen[2]);
+ // boost into lab-frame
+ TVector3 boost=sum.BoostVector();
+ // decay in rest-frame
+
+ TLorentzVector p[3];
+ double m=MASS1;
+ double costh;
+ do {
+ double r=rnd->Uniform(-1.,1.);
+ costh=r*(1.+DECAY_A*r*r);
+ } while(fabs(costh)>=1.0);
+ double phi=rnd->Uniform(-M_PI,M_PI);
+ double e=0.5*sum.M();
+ double ptot=TMath::Sqrt(e+m)*TMath::Sqrt(e-m);
+ double pz=ptot*costh;
+ double pt=TMath::Sqrt(ptot+pz)*TMath::Sqrt(ptot-pz);
+ double px=pt*cos(phi);
+ double py=pt*sin(phi);
+ p[0].SetXYZT(px,py,pz,e);
+ p[1].SetXYZT(-px,-py,-pz,e);
+ for(int i=0;i<2;i++) {
+ p[i].Boost(boost);
+ }
+ p[2]=p[0]+p[1];
+ for(int i=0;i<3;i++) {
+ fPtGen[i]=p[i].Pt();
+ fEtaGen[i]=p[i].Eta();
+ fPhiGen[i]=p[i].Phi();
+ fMGen[i]=p[i].M();
+ }
+}
+
+void ToyEvent7::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
+ for(int i=0;i<3;i++) {
+ fPtGen[i]=0.0;
+ fEtaGen[i]=0.0;
+ fPhiGen[i]=0.0;
+ }
+ TLorentzVector p[3];
+ for(int i=0;i<2;i++) {
+ p[i].SetPtEtaPhiM(rnd->Exp(15.0),rnd->Uniform(-3.,3.),
+ rnd->Uniform(-M_PI,M_PI),isData ? MASS1 : MASS1);
+ }
+ p[2]=p[0]+p[1];
+ for(int i=0;i<3;i++) {
+ fPtRec[i]=p[i].Pt();
+ fEtaRec[i]=p[i].Eta();
+ fPhiRec[i]=p[i].Phi();
+ fMRec[i]=p[i].M();
+ }
+}
+
+void ToyEvent7::GenerateReco(TRandom *rnd) {
+ if(fIsSignal) {
+ TLorentzVector p[3];
+ for(int i=0;i<2;i++) {
+ Double_t expEta=TMath::Exp(fEtaGen[i]);
+ Double_t coshEta=(expEta+1./expEta);
+ Double_t eGen=fPtGen[i]*coshEta;
+ Double_t sigmaE=
+ 0.1*TMath::Sqrt(eGen)
+ +1.0*coshEta
+ +0.01*eGen;
+ Double_t eRec;
+ do {
+ eRec=rnd->Gaus(eGen,sigmaE);
+ } while(eRec<=0.0);
+ Double_t sigmaEta=0.1+0.02*TMath::Abs(fEtaGen[i]);
+ p[i].SetPtEtaPhiM(eRec/(expEta+1./expEta),
+ rnd->Gaus(fEtaGen[i],sigmaEta),
+ remainder(rnd->Gaus(fPhiGen[i],0.03),2.*M_PI),
+ MASS1);
+ }
+ p[2]=p[0]+p[1];
+ for(int i=0;i<3;i++) {
+ fPtRec[i]=p[i].Pt();
+ fEtaRec[i]=p[i].Eta();
+ fPhiRec[i]=p[i].Phi();
+ fMRec[i]=p[i].M();
+ }
+ }
+ if(fIsSignal) {
+ do {
+ Double_t tauDiscr=0.08-0.04/(1.+fPtRec[2]/10.0);
+ Double_t sigmaDiscr=0.01;
+ fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
+ } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
+ } else {
+ do {
+ Double_t tauDiscr=0.15-0.05/(1.+fPtRec[2]/5.0)+0.1*fEtaRec[2];
+ Double_t sigmaDiscr=0.02+0.01*fEtaRec[2];
+ fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
+ } while((fDiscriminator<=0.)||(fDiscriminator>=1.));
+ }
+ fIsTriggered=false;
+ for(int i=0;i<2;i++) {
+ if(rnd->Uniform()<0.92/(TMath::Exp(-(fPtRec[i]-15.5)/2.5)+1.)) fIsTriggered=true;
+ }
+}
diff --git a/tutorials/unfold/testUnfold7b.C b/tutorials/unfold/testUnfold7b.C
new file mode 100644
index 00000000000..00ea583fc16
--- /dev/null
+++ b/tutorials/unfold/testUnfold7b.C
@@ -0,0 +1,302 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning
+///
+/// A toy test of the TUnfold package
+///
+/// This example is documented in conference proceedings:
+///
+/// arXiv:1611.01927
+/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
+///
+/// This is an example of unfolding a two-dimensional distribution
+/// also using an auxiliary measurement to constrain some background
+///
+/// The example comprises several macros
+/// - testUnfold7a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - testUnfold7b.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold7binning.xml
+/// testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - write testUnfold7_histograms.root
+///
+/// - testUnfold7c.C run the unfolding
+/// - read testUnfold7_histograms.root
+/// - write testUnfold7_result.root
+/// testUnfold7_result.ps
+///
+/// \macro_output
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+
+/* below is the content of the file testUnfold7binning.xml,
+ which is required as input to run this example.
+
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
+<TUnfoldBinning>
+<BinningNode name="fine" firstbin="1" factor="1">
+ <Axis name="pt" lowEdge="0.">
+ <Bin repeat="20" width="1."/>
+ <Bin repeat="12" width="2.5"/>
+ <Bin location="overflow" width="10"/>
+ </Axis>
+</BinningNode>
+<BinningNode name="coarse" firstbin="1" factor="1">
+ <Axis name="pt" lowEdge="0.">
+ <Bin repeat="10" width="2"/>
+ <Bin repeat="6" width="5"/>
+ <Bin location="overflow" width="10"/>
+ </Axis>
+</BinningNode>
+</TUnfoldBinning>
+
+ */
+
+#include <iostream>
+#include <map>
+#include <cmath>
+#include <TMath.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TH1.h>
+#include <TDOMParser.h>
+#include <TXMLDocument.h>
+#include "TUnfoldBinningXML.h"
+
+using namespace std;
+
+
+void testUnfold7b()
+{
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ //=======================================================
+ // Step 1: open file to save histograms and binning schemes
+
+ TFile *outputFile=new TFile("testUnfold7_histograms.root","recreate");
+
+ //=======================================================
+ // Step 2: read binning from XML
+ // and save them to output file
+
+ TUnfoldBinning *fineBinningRoot,*coarseBinningRoot;
+
+ outputFile->cd();
+
+ // read binning schemes in XML format
+
+ TDOMParser parser;
+ Int_t error=parser.ParseFile("testUnfold7binning.xml");
+ if(error) {
+ cout<<"error="<<error<<" from TDOMParser\n";
+ cout<<"==============================================================\n";
+ cout<<"Maybe the file testUnfold7binning.xml is missing?\n";
+ cout<<"The content of the file is included in the comments section\n";
+ cout<<"of this macro \"testUnfold7b.C\"\n";
+ cout<<"==============================================================\n";
+ }
+ TXMLDocument const *XMLdocument=parser.GetXMLDocument();
+ fineBinningRoot=
+ TUnfoldBinningXML::ImportXML(XMLdocument,"fine");
+ coarseBinningRoot=
+ TUnfoldBinningXML::ImportXML(XMLdocument,"coarse");
+
+ // write binning schemes to output file
+ fineBinningRoot->Write();
+ coarseBinningRoot->Write();
+
+ if(fineBinningRoot) {
+ fineBinningRoot->PrintStream(cout);
+ } else {
+ cout<<"could not read 'detector' binning\n";
+ }
+ if(coarseBinningRoot) {
+ coarseBinningRoot->PrintStream(cout);
+ } else {
+ cout<<"could not read 'generator' binning\n";
+ }
+
+ TUnfoldBinning const *fineBinning=fineBinningRoot;//->FindNode("ptfine");
+ TUnfoldBinning const *coarseBinning=coarseBinningRoot;//->FindNode("ptcoarse");
+
+
+ //=======================================================
+ // Step 3: book and fill data histograms
+
+ Float_t ptRec[3],ptGen[3],weight;
+ Int_t isTriggered,isSignal;
+
+ outputFile->cd();
+
+ TH1 *histDataRecF=fineBinning->CreateHistogram("histDataRecF");
+ TH1 *histDataRecC=coarseBinning->CreateHistogram("histDataRecC");
+ TH1 *histDataBgrF=fineBinning->CreateHistogram("histDataBgrF");
+ TH1 *histDataBgrC=coarseBinning->CreateHistogram("histDataBgrC");
+ TH1 *histDataGen=coarseBinning->CreateHistogram("histDataGen");
+
+ TFile *dataFile=new TFile("testUnfold7_data.root");
+ TTree *dataTree=(TTree *) dataFile->Get("data");
+
+ if(!dataTree) {
+ cout<<"could not read 'data' tree\n";
+ }
+
+ dataTree->ResetBranchAddresses();
+ dataTree->SetBranchAddress("ptrec",ptRec);
+ //dataTree->SetBranchAddress("discr",&discr);
+ // for real data, only the triggered events are available
+ dataTree->SetBranchAddress("istriggered",&isTriggered);
+ // data truth parameters
+ dataTree->SetBranchAddress("ptgen",ptGen);
+ dataTree->SetBranchAddress("issignal",&isSignal);
+ dataTree->SetBranchStatus("*",1);
+
+ cout<<"loop over data events\n";
+
+#define VAR_REC (ptRec[2])
+#define VAR_GEN (ptGen[2])
+
+ for(Int_t ievent=0;ievent<dataTree->GetEntriesFast();ievent++) {
+ if(dataTree->GetEntry(ievent)<=0) break;
+ // fill histogram with reconstructed quantities
+ if(isTriggered) {
+ int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
+ int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
+ histDataRecF->Fill(binF);
+ histDataRecC->Fill(binC);
+ if(!isSignal) {
+ histDataBgrF->Fill(binF);
+ histDataBgrC->Fill(binC);
+ }
+ }
+ // fill histogram with data truth parameters
+ if(isSignal) {
+ int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
+ histDataGen->Fill(binGen);
+ }
+ }
+
+ delete dataTree;
+ delete dataFile;
+
+ //=======================================================
+ // Step 4: book and fill histogram of migrations
+ // it receives events from both signal MC and background MC
+
+ outputFile->cd();
+
+ TH2 *histMcsigGenRecF=TUnfoldBinning::CreateHistogramOfMigrations
+ (coarseBinning,fineBinning,"histMcsigGenRecF");
+ TH2 *histMcsigGenRecC=TUnfoldBinning::CreateHistogramOfMigrations
+ (coarseBinning,coarseBinning,"histMcsigGenRecC");
+ TH1 *histMcsigRecF=fineBinning->CreateHistogram("histMcsigRecF");
+ TH1 *histMcsigRecC=coarseBinning->CreateHistogram("histMcsigRecC");
+ TH1 *histMcsigGen=coarseBinning->CreateHistogram("histMcsigGen");
+
+ TFile *signalFile=new TFile("testUnfold7_signal.root");
+ TTree *signalTree=(TTree *) signalFile->Get("signal");
+
+ if(!signalTree) {
+ cout<<"could not read 'signal' tree\n";
+ }
+
+ signalTree->ResetBranchAddresses();
+ signalTree->SetBranchAddress("ptrec",&ptRec);
+ //signalTree->SetBranchAddress("discr",&discr);
+ signalTree->SetBranchAddress("istriggered",&isTriggered);
+ signalTree->SetBranchAddress("ptgen",&ptGen);
+ signalTree->SetBranchAddress("weight",&weight);
+ signalTree->SetBranchStatus("*",1);
+
+ cout<<"loop over MC signal events\n";
+
+ for(Int_t ievent=0;ievent<signalTree->GetEntriesFast();ievent++) {
+ if(signalTree->GetEntry(ievent)<=0) break;
+
+ int binC=0,binF=0;
+ if(isTriggered) {
+ binF=fineBinning->GetGlobalBinNumber(VAR_REC);
+ binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
+ }
+ int binGen=coarseBinning->GetGlobalBinNumber(VAR_GEN);
+ histMcsigGenRecF->Fill(binGen,binF,weight);
+ histMcsigGenRecC->Fill(binGen,binC,weight);
+ histMcsigRecF->Fill(binF,weight);
+ histMcsigRecC->Fill(binC,weight);
+ histMcsigGen->Fill(binGen,weight);
+ }
+
+ delete signalTree;
+ delete signalFile;
+
+ outputFile->cd();
+
+ TH1 *histMcbgrRecF=fineBinning->CreateHistogram("histMcbgrRecF");
+ TH1 *histMcbgrRecC=coarseBinning->CreateHistogram("histMcbgrRecC");
+
+ TFile *bgrFile=new TFile("testUnfold7_background.root");
+ TTree *bgrTree=(TTree *) bgrFile->Get("background");
+
+ if(!bgrTree) {
+ cout<<"could not read 'background' tree\n";
+ }
+
+ bgrTree->ResetBranchAddresses();
+ bgrTree->SetBranchAddress("ptrec",&ptRec);
+ //bgrTree->SetBranchAddress("discr",&discr);
+ bgrTree->SetBranchAddress("istriggered",&isTriggered);
+ bgrTree->SetBranchAddress("weight",&weight);
+ bgrTree->SetBranchStatus("*",1);
+
+ cout<<"loop over MC background events\n";
+
+ for(Int_t ievent=0;ievent<bgrTree->GetEntriesFast();ievent++) {
+ if(bgrTree->GetEntry(ievent)<=0) break;
+
+ // here, for background only reconstructed quantities are known
+ // and only the reconstructed events are relevant
+ if(isTriggered) {
+ int binF=fineBinning->GetGlobalBinNumber(VAR_REC);
+ int binC=coarseBinning->GetGlobalBinNumber(VAR_REC);
+ histMcbgrRecF->Fill(binF,weight);
+ histMcbgrRecC->Fill(binC,weight);
+ }
+ }
+
+ delete bgrTree;
+ delete bgrFile;
+
+ outputFile->Write();
+ delete outputFile;
+
+}
diff --git a/tutorials/unfold/testUnfold7binning.xml b/tutorials/unfold/testUnfold7binning.xml
new file mode 100644
index 00000000000..635179045b8
--- /dev/null
+++ b/tutorials/unfold/testUnfold7binning.xml
@@ -0,0 +1,18 @@
+<?xml version="1.0" encoding="UTF-8" standalone="no"?>
+<!DOCTYPE TUnfoldBinning SYSTEM "tunfoldbinning.dtd">
+<TUnfoldBinning>
+<BinningNode name="fine" firstbin="1" factor="1">
+ <Axis name="pt" lowEdge="0.">
+ <Bin repeat="20" width="1."/>
+ <Bin repeat="12" width="2.5"/>
+ <Bin location="overflow" width="10"/>
+ </Axis>
+</BinningNode>
+<BinningNode name="coarse" firstbin="1" factor="1">
+ <Axis name="pt" lowEdge="0.">
+ <Bin repeat="10" width="2"/>
+ <Bin repeat="6" width="5"/>
+ <Bin location="overflow" width="10"/>
+ </Axis>
+</BinningNode>
+</TUnfoldBinning>
diff --git a/tutorials/unfold/testUnfold7c.C b/tutorials/unfold/testUnfold7c.C
new file mode 100644
index 00000000000..e889e99f324
--- /dev/null
+++ b/tutorials/unfold/testUnfold7c.C
@@ -0,0 +1,1655 @@
+/// \file
+/// \ingroup tutorial_unfold
+/// \notebook
+/// Test program for the classes TUnfoldDensity and TUnfoldBinning.
+///
+/// A toy test of the TUnfold package
+///
+///
+/// This example is documented in conference proceedings:
+///
+/// arXiv:1611.01927
+/// 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
+///
+/// This is an example of unfolding a one-dimensional distribution. It compares
+/// various unfolding methods:
+///
+/// matrix inversion, template fit, L-curve scan,
+/// iterative unfolding, etc
+///
+/// Further details can be found in talk by S.Schmitt at:
+///
+/// XII Quark Confinement and the Hadron Spectrum
+/// 29.8. - 3.9.2016 Thessaloniki, Greece
+/// statictics session (+proceedings)
+///
+/// The example comprises several macros
+/// - testUnfold7a.C create root files with TTree objects for
+/// signal, background and data
+/// - write files testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - testUnfold7b.C loop over trees and fill histograms based on the
+/// TUnfoldBinning objects
+/// - read testUnfold7binning.xml
+/// testUnfold7_signal.root
+/// testUnfold7_background.root
+/// testUnfold7_data.root
+///
+/// - write testUnfold7_histograms.root
+///
+/// - testUnfold7c.C run the unfolding
+/// - read testUnfold7_histograms.root
+/// - write testUnfold7_result.root
+/// - write many histograms, to compare various unfolding methods
+///
+/// \macro_output
+/// \macro_image
+/// \macro_code
+///
+/// **Version 17.6, in parallel to changes in TUnfold**
+///
+/// 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/>.
+///
+/// \author Stefan Schmitt DESY, 14.10.2008
+
+#include <iostream>
+#include <cmath>
+#include <map>
+#include <TMath.h>
+#include <TCanvas.h>
+#include <TStyle.h>
+#include <TGraph.h>
+#include <TGraphErrors.h>
+#include <TFile.h>
+#include <TROOT.h>
+#include <TText.h>
+#include <TLine.h>
+#include <TLegend.h>
+#include <TH1.h>
+#include <TF1.h>
+#include <TFitter.h>
+#include <TMatrixD.h>
+#include <TMatrixDSym.h>
+#include <TVectorD.h>
+#include <TMatrixDSymEigen.h>
+#include <TFitResult.h>
+#include <TRandom3.h>
+#include "TUnfoldDensity.h"
+
+using namespace std;
+
+// #define PRINT_MATRIX_L
+
+#define TEST_INPUT_COVARIANCE
+
+void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning);
+void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX);
+
+TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY);
+TH1 *AddOverflowX(TH1 *h,double width);
+
+void DrawOverflowX(TH1 *h,double posy);
+void DrawOverflowY(TH1 *h,double posx);
+
+
+double const kLegendFontSize=0.05;
+int kNbinC=0;
+
+void DrawPadProbability(TH2 *h);
+void DrawPadEfficiency(TH1 *h);
+void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
+ TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij);
+void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
+ char const *text=0,double tau=0.0,vector<double> const *r=0,
+ TF1 *f=0);
+void DrawPadCorrelations(TH2 *h,
+ vector<pair<TF1*,vector<double> > > const *table);
+
+TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,char const *text,
+ vector<pair<TF1*,vector<double> > > &table,int niter=0);
+
+void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
+ vector<double> > > const &table,int color,
+ TGraph *graph[4],int style);
+void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
+ TH1 *hist[4],int color,int style,int fillStyle);
+
+#ifdef WITH_IDS
+void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr);
+
+void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_,TMatrixD *Am_,
+ Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,
+ TVectorD* &unfres2IDS_ ,TVectorD *&soustr);
+#endif
+
+TRandom3 *g_rnd;
+
+void testUnfold7c()
+{
+ gErrorIgnoreLevel=kInfo;
+ // switch on histogram errors
+ TH1::SetDefaultSumw2();
+
+ gStyle->SetOptStat(0);
+
+ g_rnd=new TRandom3(4711);
+
+ //==============================================
+ // step 1 : open output file
+ TFile *outputFile=new TFile("testUnfold7_results.root","recreate");
+
+ //==============================================
+ // step 2 : read binning schemes and input histograms
+ TFile *inputFile=new TFile("testUnfold7_histograms.root");
+
+ outputFile->cd();
+
+ TUnfoldBinning *fineBinning,*coarseBinning;
+
+ inputFile->GetObject("fine",fineBinning);
+ inputFile->GetObject("coarse",coarseBinning);
+
+ if((!fineBinning)||(!coarseBinning)) {
+ cout<<"problem to read binning schemes\n";
+ }
+
+ // save binning schemes to output file
+ fineBinning->Write();
+ coarseBinning->Write();
+
+ // read histograms
+#define READ(TYPE,binning,name) \
+ TYPE *name[3]; inputFile->GetObject(#name,name[0]); \
+ name[0]->Write(); \
+ if(!name[0]) cout<<"Error reading " #name "\n"; \
+ CreateHistogramCopies(name,binning);
+
+ outputFile->cd();
+
+ READ(TH1,fineBinning,histDataRecF);
+ READ(TH1,coarseBinning,histDataRecC);
+ READ(TH1,fineBinning,histDataBgrF);
+ READ(TH1,coarseBinning,histDataBgrC);
+ READ(TH1,coarseBinning,histDataGen);
+
+ READ(TH2,fineBinning,histMcsigGenRecF);
+ READ(TH2,coarseBinning,histMcsigGenRecC);
+ READ(TH1,fineBinning,histMcsigRecF);
+ READ(TH1,coarseBinning,histMcsigRecC);
+ READ(TH1,coarseBinning,histMcsigGen);
+
+ READ(TH1,fineBinning,histMcbgrRecF);
+ READ(TH1,coarseBinning,histMcbgrRecC);
+
+ TH1 *histOutputCtau0[3];
+ TH2 *histRhoCtau0;
+ TH1 *histOutputCLCurve[3];
+ //TH2 *histRhoCLCurve;
+ TH2 *histProbC;
+ double tauMin=1.e-4;
+ double tauMax=1.e-1;
+ double fBgr=1.0; // 0.2/0.25;
+ double biasScale=1.0;
+ TUnfold::ERegMode mode=TUnfold::kRegModeSize; //Derivative;
+
+ //double tauC;
+ {
+ TUnfoldDensity *tunfoldC=
+ new TUnfoldDensity(histMcsigGenRecC[0],
+ TUnfold::kHistMapOutputHoriz,
+ mode,
+ TUnfold::kEConstraintNone,//Area,
+ TUnfoldDensity::kDensityModeNone,
+ coarseBinning,
+ coarseBinning);
+ tunfoldC->SetInput(histDataRecC[0],biasScale);
+ tunfoldC->SubtractBackground(histMcbgrRecC[0],"BGR",fBgr,0.0);
+ tunfoldC->DoUnfold(0.);
+ histOutputCtau0[0]=tunfoldC->GetOutput("histOutputCtau0");
+ histRhoCtau0=tunfoldC->GetRhoIJtotal("histRhoCtau0");
+ CreateHistogramCopies(histOutputCtau0,coarseBinning);
+ tunfoldC->ScanLcurve(50,tauMin,tauMax,0);
+ /* tauC= */tunfoldC->GetTau();
+ //tunfoldC->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
+ histOutputCLCurve[0]=tunfoldC->GetOutput("histOutputCLCurve");
+ /* histRhoCLCurve= */tunfoldC->GetRhoIJtotal("histRhoCLCurve");
+ CreateHistogramCopies(histOutputCLCurve,coarseBinning);
+ histProbC=tunfoldC->GetProbabilityMatrix("histProbC",";P_T(gen);P_T(rec)");
+ }
+ TH1 *histOutputFtau0[3];
+ TH2 *histRhoFtau0;
+ TH1 *histOutputFLCurve[3];
+ //TH2 *histRhoFLCurve;
+ TH2 *histProbF;
+ TGraph *lCurve;
+ TSpline *logTauX,*logTauY;
+ tauMin=3.E-4;
+ tauMax=3.E-2;
+ //double tauF;
+ {
+ TUnfoldDensity *tunfoldF=
+ new TUnfoldDensity(histMcsigGenRecF[0],
+ TUnfold::kHistMapOutputHoriz,
+ mode,
+ TUnfold::kEConstraintNone,//Area,
+ TUnfoldDensity::kDensityModeNone,
+ coarseBinning,
+ fineBinning);
+ tunfoldF->SetInput(histDataRecF[0],biasScale);
+ tunfoldF->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
+ tunfoldF->DoUnfold(0.);
+ histOutputFtau0[0]=tunfoldF->GetOutput("histOutputFtau0");
+ histRhoFtau0=tunfoldF->GetRhoIJtotal("histRhoFtau0");
+ CreateHistogramCopies(histOutputFtau0,coarseBinning);
+ tunfoldF->ScanLcurve(50,tauMin,tauMax,0);
+ //tunfoldF->DoUnfold(tauC);
+ /* tauF= */tunfoldF->GetTau();
+ //tunfoldF->ScanTau(50,1.E-7,1.E-1,0,TUnfoldDensity::kEScanTauRhoAvg);
+ histOutputFLCurve[0]=tunfoldF->GetOutput("histOutputFLCurve");
+ /* histRhoFLCurve= */tunfoldF->GetRhoIJtotal("histRhoFLCurve");
+ CreateHistogramCopies(histOutputFLCurve,coarseBinning);
+ histProbF=tunfoldF->GetProbabilityMatrix("histProbF",";P_T(gen);P_T(rec)");
+ }
+ TH1 *histOutputFAtau0[3];
+ TH2 *histRhoFAtau0;
+ TH1 *histOutputFALCurve[3];
+ TH2 *histRhoFALCurve;
+ TH1 *histOutputFArho[3];
+ TH2 *histRhoFArho;
+ TSpline *rhoScan=0;
+ TSpline *logTauCurvature=0;
+
+ double tauFA,tauFArho;
+ {
+ TUnfoldDensity *tunfoldFA=
+ new TUnfoldDensity(histMcsigGenRecF[0],
+ TUnfold::kHistMapOutputHoriz,
+ mode,
+ TUnfold::kEConstraintArea,
+ TUnfoldDensity::kDensityModeNone,
+ coarseBinning,
+ fineBinning);
+ tunfoldFA->SetInput(histDataRecF[0],biasScale);
+ tunfoldFA->SubtractBackground(histMcbgrRecF[0],"BGR",fBgr,0.0);
+ tunfoldFA->DoUnfold(0.);
+ histOutputFAtau0[0]=tunfoldFA->GetOutput("histOutputFAtau0");
+ histRhoFAtau0=tunfoldFA->GetRhoIJtotal("histRhoFAtau0");
+ CreateHistogramCopies(histOutputFAtau0,coarseBinning);
+ tunfoldFA->ScanTau(50,tauMin,tauMax,&rhoScan,TUnfoldDensity::kEScanTauRhoAvg);
+ tauFArho=tunfoldFA->GetTau();
+ histOutputFArho[0]=tunfoldFA->GetOutput("histOutputFArho");
+ histRhoFArho=tunfoldFA->GetRhoIJtotal("histRhoFArho");
+ CreateHistogramCopies(histOutputFArho,coarseBinning);
+
+ tunfoldFA->ScanLcurve(50,tauMin,tauMax,&lCurve,&logTauX,&logTauY,&logTauCurvature);
+ tauFA=tunfoldFA->GetTau();
+ histOutputFALCurve[0]=tunfoldFA->GetOutput("histOutputFALCurve");
+ histRhoFALCurve=tunfoldFA->GetRhoIJtotal("histRhoFALCurve");
+ CreateHistogramCopies(histOutputFALCurve,coarseBinning);
+ }
+ lCurve->Write();
+ logTauX->Write();
+ logTauY->Write();
+
+
+ double widthC=coarseBinning->GetBinSize(histProbC->GetNbinsY()+1);
+ double widthF=fineBinning->GetBinSize(histProbF->GetNbinsY()+1);
+
+ TH2 *histProbCO=AddOverflowXY(histProbC,widthC,widthC);
+ TH2 *histProbFO=AddOverflowXY(histProbF,widthC,widthF);
+
+ // efficiency
+ TH1 *histEfficiencyC=histProbCO->ProjectionX("histEfficiencyC");
+ kNbinC=histProbCO->GetNbinsX();
+
+ // reconstructed quantities with overflow (coarse binning)
+ // MC: add signal and bgr
+ TH1 *histMcsigRecCO=AddOverflowX(histMcsigRecC[2],widthC);
+ TH1 *histMcbgrRecCO=AddOverflowX(histMcbgrRecC[2],widthC);
+ histMcbgrRecCO->Scale(fBgr);
+ TH1 *histMcRecCO=(TH1 *)histMcsigRecCO->Clone("histMcRecC0");
+ histMcRecCO->Add(histMcsigRecCO,histMcbgrRecCO);
+ TH1 *histDataRecCO=AddOverflowX(histDataRecC[2],widthC);
+ //TH1 *histDataRecCNO=AddOverflowX(histDataRecC[1],widthC);
+
+ TH1 *histMcsigRecFO=AddOverflowX(histMcsigRecF[2],widthF);
+ TH1 *histMcbgrRecFO=AddOverflowX(histMcbgrRecF[2],widthF);
+ histMcbgrRecFO->Scale(fBgr);
+ TH1 *histMcRecFO=(TH1 *)histMcsigRecFO->Clone("histMcRecF0");
+ histMcRecFO->Add(histMcsigRecFO,histMcbgrRecFO);
+ TH1 *histDataRecFO=AddOverflowX(histDataRecF[2],widthF);
+
+ // truth level with overflow
+ TH1 *histMcsigGenO=AddOverflowX(histMcsigGen[2],widthC);
+ TH1 *histDataGenO=AddOverflowX(histDataGen[2],widthC);
+
+ // unfolding result with overflow
+ TH1 *histOutputCtau0O=AddOverflowX(histOutputCtau0[2],widthC);
+ TH2 *histRhoCtau0O=AddOverflowXY(histRhoCtau0,widthC,widthC);
+ //TH1 *histOutputCLCurveO=AddOverflowX(histOutputCLCurve[2],widthC);
+ //TH2 *histRhoCLCurveO=AddOverflowXY(histRhoCLCurve,widthC,widthC);
+ TH1 *histOutputFtau0O=AddOverflowX(histOutputFtau0[2],widthC);
+ TH2 *histRhoFtau0O=AddOverflowXY(histRhoFtau0,widthC,widthC);
+ TH1 *histOutputFAtau0O=AddOverflowX(histOutputFAtau0[2],widthC);
+ TH2 *histRhoFAtau0O=AddOverflowXY(histRhoFAtau0,widthC,widthC);
+ //TH1 *histOutputFLCurveO=AddOverflowX(histOutputFLCurve[2],widthC);
+ //TH2 *histRhoFLCurveO=AddOverflowXY(histRhoFLCurve,widthC,widthC);
+ TH1 *histOutputFALCurveO=AddOverflowX(histOutputFALCurve[2],widthC);
+ TH2 *histRhoFALCurveO=AddOverflowXY(histRhoFALCurve,widthC,widthC);
+ TH1 *histOutputFArhoO=AddOverflowX(histOutputFArho[2],widthC);
+ TH2 *histRhoFArhoO=AddOverflowXY(histRhoFArho,widthC,widthC);
+
+ // bin-by-bin
+ TH2 *histRhoBBBO=(TH2 *)histRhoCtau0O->Clone("histRhoBBBO");
+ for(int i=1;i<=histRhoBBBO->GetNbinsX();i++) {
+ for(int j=1;j<=histRhoBBBO->GetNbinsX();j++) {
+ histRhoBBBO->SetBinContent(i,j,(i==j)?1.:0.);
+ }
+ }
+ TH1 *histDataBgrsub=(TH1 *)histDataRecCO->Clone("histDataBgrsub");
+ histDataBgrsub->Add(histMcbgrRecCO,-fBgr);
+ TH1 *histOutputBBBO=(TH1 *)histDataBgrsub->Clone("histOutputBBBO");
+ histOutputBBBO->Divide(histMcsigRecCO);
+ histOutputBBBO->Multiply(histMcsigGenO);
+
+ // iterative
+ int niter=1000;
+ cout<<"maximum number of iterations: "<<niter<<"\n";
+
+ vector <TH1 *>histOutputAgo,histOutputAgorep;
+ vector <TH2 *>histRhoAgo,histRhoAgorep;
+ vector<int> nIter;
+ histOutputAgo.push_back((TH1*)histMcsigGenO->Clone("histOutputAgo-1"));
+ histOutputAgorep.push_back((TH1*)histMcsigGenO->Clone("histOutputAgorep-1"));
+ histRhoAgo.push_back((TH2*)histRhoBBBO->Clone("histRhoAgo-1"));
+ histRhoAgorep.push_back((TH2*)histRhoBBBO->Clone("histRhoAgorep-1"));
+ nIter.push_back(-1);
+
+ int nx=histProbCO->GetNbinsX();
+ int ny=histProbCO->GetNbinsY();
+ TMatrixD covAgo(nx+ny,nx+ny);
+ TMatrixD A(ny,nx);
+ TMatrixD AToverEps(nx,ny);
+ for(int i=0;i<nx;i++) {
+ double epsilonI=0.;
+ for(int j=0;j<ny;j++) {
+ epsilonI+= histProbCO->GetBinContent(i+1,j+1);
+ }
+ for(int j=0;j<ny;j++) {
+ double aji=histProbCO->GetBinContent(i+1,j+1);
+ A(j,i)=aji;
+ AToverEps(i,j)=aji/epsilonI;
+ }
+ }
+ for(int i=0;i<nx;i++) {
+ covAgo(i,i)=TMath::Power
+ (histOutputAgo[0]->GetBinError(i+1)
+ *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1),2.);
+ }
+ for(int i=0;i<ny;i++) {
+ covAgo(i+nx,i+nx)=TMath::Power
+ (histDataRecCO->GetBinError(i+1)
+ *histDataRecCO->GetXaxis()->GetBinWidth(i+1),2.);
+ }
+#define NREPLICA 300
+ vector<TVectorD *> y(NREPLICA);
+ vector<TVectorD *> yMb(NREPLICA);
+ vector<TVectorD *> yErr(NREPLICA);
+ vector<TVectorD *> x(NREPLICA);
+ TVectorD b(nx);
+ for(int nr=0;nr<NREPLICA;nr++) {
+ x[nr]=new TVectorD(nx);
+ y[nr]=new TVectorD(ny);
+ yMb[nr]=new TVectorD(ny);
+ yErr[nr]=new TVectorD(ny);
+ }
+ for(int i=0;i<nx;i++) {
+ (*x[0])(i)=histOutputAgo[0]->GetBinContent(i+1)
+ *histOutputAgo[0]->GetXaxis()->GetBinWidth(i+1);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ (*x[nr])(i)=(*x[0])(i);
+ }
+ }
+ for(int i=0;i<ny;i++) {
+ (*y[0])(i)=histDataRecCO->GetBinContent(i+1)
+ *histDataRecCO->GetXaxis()->GetBinWidth(i+1);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ (*y[nr])(i)=g_rnd->Poisson((*y[0])(i));
+ (*yErr[nr])(i)=TMath::Sqrt((*y[nr])(i));
+ }
+ b(i)=histMcbgrRecCO->GetBinContent(i+1)*
+ histMcbgrRecCO->GetXaxis()->GetBinWidth(i+1);
+ for(int nr=0;nr<NREPLICA;nr++) {
+ (*yMb[nr])(i)=(*y[nr])(i)-b(i);
+ }
+ }
+ for(int iter=0;iter<=niter;iter++) {
+ if(!(iter %100)) cout<<iter<<"\n";
+ for(int nr=0;nr<NREPLICA;nr++) {
+ TVectorD yrec=A*(*x[nr])+b;
+ TVectorD yOverYrec(ny);
+ for(int j=0;j<ny;j++) {
+ yOverYrec(j)=(*y[nr])(j)/yrec(j);
+ }
+ TVectorD f=AToverEps * yOverYrec;
+ TVectorD xx(nx);
+ for(int i=0;i<nx;i++) {
+ xx(i) = (*x[nr])(i) * f(i);
+ }
+ if(nr==0) {
+ TMatrixD xdf_dr=AToverEps;
+ for(int i=0;i<nx;i++) {
+ for(int j=0;j<ny;j++) {
+ xdf_dr(i,j) *= (*x[nr])(i);
+ }
+ }
+ TMatrixD dr_dxdy(ny,nx+ny);
+ for(int j=0;j<ny;j++) {
+ dr_dxdy(j,nx+j)=1.0/yrec(j);
+ for(int i=0;i<nx;i++) {
+ dr_dxdy(j,i)= -yOverYrec(j)/yrec(j)*A(j,i);
+ }
+ }
+ TMatrixD dxy_dxy(nx+ny,nx+ny);
+ dxy_dxy.SetSub(0,0,xdf_dr*dr_dxdy);
+ for(int i=0;i<nx;i++) {
+ dxy_dxy(i,i) +=f(i);
+ }
+ for(int i=0;i<ny;i++) {
+ dxy_dxy(nx+i,nx+i) +=1.0;
+ }
+ TMatrixD VDT(covAgo,TMatrixD::kMultTranspose,dxy_dxy);
+ covAgo= dxy_dxy*VDT;
+ }
+ (*x[nr])=xx;
+ }
+ if((iter<=25)||
+ ((iter<=100)&&(iter %5==0))||
+ ((iter<=1000)&&(iter %50==0))||
+ (iter %1000==0)) {
+ nIter.push_back(iter);
+ TH1 * h=(TH1*)histOutputAgo[0]->Clone
+ (TString::Format("histOutputAgo%d",iter));
+ histOutputAgo.push_back(h);
+ for(int i=0;i<nx;i++) {
+ double bw=h->GetXaxis()->GetBinWidth(i+1);
+ h->SetBinContent(i+1,(*x[0])(i)/bw);
+ h->SetBinError(i+1,TMath::Sqrt(covAgo(i,i))/bw);
+ }
+ TH2 *h2=(TH2*)histRhoAgo[0]->Clone
+ (TString::Format("histRhoAgo%d",iter));
+ histRhoAgo.push_back(h2);
+ for(int i=0;i<nx;i++) {
+ for(int j=0;j<nx;j++) {
+ double rho= covAgo(i,j)/TMath::Sqrt(covAgo(i,i)*covAgo(j,j));
+ if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
+ cout<<"bad error matrix: iter="<<iter<<"\n";
+ exit(0);
+ }
+ h2->SetBinContent(i+1,j+1,rho);
+ }
+ }
+ // error and correlations from replica analysis
+ h=(TH1*)histOutputAgo[0]->Clone
+ (TString::Format("histOutputAgorep%d",iter));
+ h2=(TH2*)histRhoAgo[0]->Clone
+ (TString::Format("histRhoAgorep%d",iter));
+ histOutputAgorep.push_back(h);
+ histRhoAgorep.push_back(h2);
+
+ TVectorD mean(nx);
+ double w=1./(NREPLICA-1.);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ mean += w* *x[nr];
+ }
+ TMatrixD covAgorep(nx,nx);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ //TMatrixD dx= (*x)-mean;
+ TMatrixD dx(nx,1);
+ for(int i=0;i<nx;i++) {
+ dx(i,0)= (*x[nr])(i)-(*x[0])(i);
+ }
+ covAgorep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
+ }
+
+ for(int i=0;i<nx;i++) {
+ double bw=h->GetXaxis()->GetBinWidth(i+1);
+ h->SetBinContent(i+1,(*x[0])(i)/bw);
+ h->SetBinError(i+1,TMath::Sqrt(covAgorep(i,i))/bw);
+ // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
+ }
+ for(int i=0;i<nx;i++) {
+ for(int j=0;j<nx;j++) {
+ double rho= covAgorep(i,j)/
+ TMath::Sqrt(covAgorep(i,i)*covAgorep(j,j));
+ if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
+ cout<<"bad error matrix: iter="<<iter<<"\n";
+ exit(0);
+ }
+ h2->SetBinContent(i+1,j+1,rho);
+ }
+ }
+ }
+ }
+
+#ifdef WITH_IDS
+ // IDS Malaescu
+ int niterIDS=100;
+ vector<TVectorD*> unfresIDS(NREPLICA),soustr(NREPLICA);
+ cout<<"IDS number of iterations: "<<niterIDS<<"\n";
+ TMatrixD *Am_IDS[NREPLICA];
+ TMatrixD A_IDS(ny,nx);
+ for(int nr=0;nr<NREPLICA;nr++) {
+ Am_IDS[nr]=new TMatrixD(ny,nx);
+ }
+ for(int iy=0;iy<ny;iy++) {
+ for(int ix=0;ix<nx;ix++) {
+ A_IDS(iy,ix)=histMcsigGenRecC[0]->GetBinContent(ix+1,iy+1);
+ }
+ }
+ double lambdaL=0.;
+ Double_t lambdaUmin = 1.0000002;
+ Double_t lambdaMmin = 0.0000001;
+ Double_t lambdaS = 0.000001;
+ double lambdaU=lambdaUmin;
+ double lambdaM=lambdaMmin;
+ vector<TH1 *> histOutputIDS;
+ vector<TH2 *> histRhoIDS;
+ histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS-1"));
+ histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS-1"));
+ histOutputIDS.push_back((TH1*)histOutputAgo[0]->Clone("histOutputIDS0"));
+ histRhoIDS.push_back((TH2*)histRhoAgo[0]->Clone("histRhoIDS0"));
+ for(int iter=1;iter<=niterIDS;iter++) {
+ if(!(iter %10)) cout<<iter<<"\n";
+
+ for(int nr=0;nr<NREPLICA;nr++) {
+ if(iter==1) {
+ IDSfirst(yMb[nr],yErr[nr],&A_IDS,lambdaL,unfresIDS[nr],soustr[nr]);
+ } else {
+ IDSiterate(yMb[nr],yErr[nr],&A_IDS,Am_IDS[nr],
+ lambdaU,lambdaM,lambdaS,
+ unfresIDS[nr],soustr[nr]);
+ }
+ }
+ unsigned ix;
+ for(ix=0;ix<nIter.size();ix++) {
+ if(nIter[ix]==iter) break;
+ }
+ if(ix<nIter.size()) {
+ TH1 * h=(TH1*)histOutputIDS[0]->Clone
+ (TString::Format("histOutputIDS%d",iter));
+ TH2 *h2=(TH2*)histRhoIDS[0]->Clone
+ (TString::Format("histRhoIDS%d",iter));
+ histOutputIDS.push_back(h);
+ histRhoIDS.push_back(h2);
+ TVectorD mean(nx);
+ double w=1./(NREPLICA-1.);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ mean += w* (*unfresIDS[nr]);
+ }
+ TMatrixD covIDSrep(nx,nx);
+ for(int nr=1;nr<NREPLICA;nr++) {
+ //TMatrixD dx= (*x)-mean;
+ TMatrixD dx(nx,1);
+ for(int i=0;i<nx;i++) {
+ dx(i,0)= (*unfresIDS[nr])(i)-(*unfresIDS[0])(i);
+ }
+ covIDSrep += w*TMatrixD(dx,TMatrixD::kMultTranspose,dx);
+ }
+ for(int i=0;i<nx;i++) {
+ double bw=h->GetXaxis()->GetBinWidth(i+1);
+ h->SetBinContent(i+1,(*unfresIDS[0])(i)/bw/
+ histEfficiencyC->GetBinContent(i+1));
+ h->SetBinError(i+1,TMath::Sqrt(covIDSrep(i,i))/bw/
+ histEfficiencyC->GetBinContent(i+1));
+ // cout<<i<<" "<<(*x[0])(i)/bw<<" +/-"<<TMath::Sqrt(covAgorep(i,i))/bw<<" "<<TMath::Sqrt(covAgo(i,i))/bw<<"\n";
+ }
+ for(int i=0;i<nx;i++) {
+ for(int j=0;j<nx;j++) {
+ double rho= covIDSrep(i,j)/
+ TMath::Sqrt(covIDSrep(i,i)*covIDSrep(j,j));
+ if((i!=j)&&(TMath::Abs(rho)>=1.0)) {
+ cout<<"bad error matrix: iter="<<iter<<"\n";
+ exit(0);
+ }
+ h2->SetBinContent(i+1,j+1,rho);
+ }
+ }
+ }
+ }
+#endif
+
+ //double NEdSmc=histDataBgrsub->GetSumOfWeights();
+
+ vector<pair<TF1 *,vector<double> > > table;
+
+ TCanvas *c1=new TCanvas("c1","",600,600);
+ TCanvas *c2sq=new TCanvas("c2sq","",600,600);
+ c2sq->Divide(1,2);
+ TCanvas *c2w=new TCanvas("c2w","",600,300);
+ c2w->Divide(2,1);
+ TCanvas *c4=new TCanvas("c4","",600,600);
+ c4->Divide(2,2);
+ //TCanvas *c3n=new TCanvas("c3n","",600,600);
+ TPad *subn[3];
+ //gROOT->SetStyle("xTimes2");
+ subn[0]= new TPad("subn0","",0.,0.5,1.,1.);
+ //gROOT->SetStyle("square");
+ subn[1]= new TPad("subn1","",0.,0.,0.5,0.5);
+ subn[2]= new TPad("subn2","",0.5,0.0,1.,0.5);
+ for(int i=0;i<3;i++) {
+ subn[i]->SetFillStyle(0);
+ subn[i]->Draw();
+ }
+ TCanvas *c3c=new TCanvas("c3c","",600,600);
+ TPad *subc[3];
+ //gROOT->SetStyle("xTimes2");
+ subc[0]= new TPad("sub0","",0.,0.5,1.,1.);
+ //gROOT->SetStyle("squareCOLZ");
+ subc[1]= new TPad("sub1","",0.,0.,0.5,0.5);
+ //gROOT->SetStyle("square");
+ subc[2]= new TPad("sub2","",0.5,0.0,1.,0.5);
+ for(int i=0;i<3;i++) {
+ subc[i]->SetFillStyle(0);
+ subc[i]->Draw();
+ }
+
+ //=========================== example ==================================
+
+ c2w->cd(1);
+ DrawPadTruth(histMcsigGenO,histDataGenO,0);
+ c2w->cd(2);
+ DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,0,0,0);
+ c2w->SaveAs("exampleTR.eps");
+
+ //=========================== example ==================================
+
+ c2w->cd(1);
+ DrawPadProbability(histProbCO);
+ c2w->cd(2);
+ DrawPadEfficiency(histEfficiencyC);
+ c2w->SaveAs("exampleAE.eps");
+
+ int iFitInversion=table.size();
+ DoFit(histOutputCtau0O,histRhoCtau0O,histDataGenO,"inversion",table);
+ //=========================== inversion ==================================
+
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputCtau0O,"inversion",0.,
+ &table[table.size()-1].second);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoCtau0O,&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
+ histOutputCtau0O,histProbCO,histRhoCtau0O);
+ c3c->SaveAs("inversion.eps");
+
+
+ DoFit(histOutputFtau0O,histRhoFtau0O,histDataGenO,"template",table);
+ //=========================== template ==================================
+
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFtau0O,"fit",0.,
+ &table[table.size()-1].second);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoFtau0O,&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
+ histOutputFtau0O,histProbFO,histRhoFtau0O);
+ c3c->SaveAs("template.eps");
+
+ DoFit(histOutputFAtau0O,histRhoFAtau0O,histDataGenO,"template+area",table);
+ //=========================== template+area ==================================
+
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFAtau0O,"fit",0.,
+ &table[table.size()-1].second);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoFAtau0O,&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
+ histOutputFAtau0O,histProbFO,histRhoFAtau0O);
+ c3c->SaveAs("templateA.eps");
+
+ int iFitFALCurve=table.size();
+ DoFit(histOutputFALCurveO,histRhoFALCurveO,histDataGenO,"Tikhonov+area",table);
+ //=========================== template+area+tikhonov =====================
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
+ &table[table.size()-1].second);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoFALCurveO,&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
+ histOutputFALCurveO,histProbFO,histRhoFALCurveO);
+ c3c->SaveAs("lcurveFA.eps");
+
+ int iFitFArho=table.size();
+ DoFit(histOutputFArhoO,histRhoFArhoO,histDataGenO,"min(rhomax)",table);
+ //=========================== template+area+tikhonov =====================
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
+ &table[table.size()-1].second);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoFArho,&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecFO,histMcbgrRecFO,histDataRecFO,
+ histOutputFArhoO,histProbFO,histRhoFArhoO);
+ c3c->SaveAs("rhoscanFA.eps");
+
+ int iFitBinByBin=table.size();
+ DoFit(histOutputBBBO,histRhoBBBO,histDataGenO,"bin-by-bin",table);
+ //=========================== bin-by-bin =================================
+ //c->cd(1);
+ //DrawPadProbability(histProbCO);
+ //c->cd(2);
+ //DrawPadCorrelations(histRhoBBBO,&table);
+ c2sq->cd(1);
+ DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
+ histOutputBBBO,histProbCO,histRhoBBBO);
+ c2sq->cd(2);
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputBBBO,"bin-by-bin",0.,
+ &table[table.size()-1].second);
+ c2sq->SaveAs("binbybin.eps");
+
+
+ //=========================== iterative ===================================
+ int iAgoFirstFit=table.size();
+ for(size_t i=1;i<histRhoAgorep.size();i++) {
+ int n=nIter[i];
+ bool isFitted=false;
+ DoFit(histOutputAgorep[i],histRhoAgorep[i],histDataGenO,
+ TString::Format("iterative, N=%d",n),table,n);
+ isFitted=true;
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputAgorep[i],
+ TString::Format("iterative N=%d",nIter[i]),0.,
+ isFitted ? &table[table.size()-1].second : 0);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoAgorep[i],&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
+ histOutputAgorep[i],histProbCO,histRhoAgorep[i]);
+ c3c->SaveAs(TString::Format("iterative%d.eps",nIter[i]));
+ }
+ int iAgoLastFit=table.size();
+
+#ifdef WITH_IDS
+ int iIDSFirstFit=table.size();
+
+ //=========================== IDS ===================================
+
+ for(size_t i=2;i<histRhoIDS.size();i++) {
+ int n=nIter[i];
+ bool isFitted=false;
+ DoFit(histOutputIDS[i],histRhoIDS[i],histDataGenO,
+ TString::Format("IDS, N=%d",n),table,n);
+ isFitted=true;
+ subc[0]->cd();
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputIDS[i],
+ TString::Format("IDS N=%d",nIter[i]),0.,
+ isFitted ? &table[table.size()-1].second : 0);
+ subc[1]->cd();
+ DrawPadCorrelations(histRhoIDS[i],&table);
+ subc[2]->cd();
+ DrawPadReco(histMcRecCO,histMcbgrRecCO,histDataRecCO,
+ histOutputIDS[i],histProbCO,histRhoIDS[i]);
+ c3c->SaveAs(TString::Format("ids%d.eps",nIter[i]));
+ }
+ int iIDSLastFit=table.size();
+#endif
+
+ int nfit=table.size();
+ TH1D *fitChindf=new TH1D("fitChindf",";algorithm;#chi^{2}/NDF",nfit,0,nfit);
+ TH1D *fitNorm=new TH1D("fitNorm",";algorithm;Landau amplitude [1/GeV]",nfit,0,nfit);
+ TH1D *fitMu=new TH1D("fitMu",";algorithm;Landau #mu [GeV]",nfit,0,nfit);
+ TH1D *fitSigma=new TH1D("fitSigma",";algorithm;Landau #sigma [GeV]",nfit,0,nfit);
+ for(int fit=0;fit<nfit;fit++) {
+ TF1 *f=table[fit].first;
+ vector<double> const &r=table[fit].second;
+ fitChindf->GetXaxis()->SetBinLabel(fit+1,f->GetName());
+ fitNorm->GetXaxis()->SetBinLabel(fit+1,f->GetName());
+ fitMu->GetXaxis()->SetBinLabel(fit+1,f->GetName());
+ fitSigma->GetXaxis()->SetBinLabel(fit+1,f->GetName());
+ double chi2=r[0];
+ double ndf=r[1];
+ fitChindf->SetBinContent(fit+1,chi2/ndf);
+ fitChindf->SetBinError(fit+1,TMath::Sqrt(2./ndf));
+ fitNorm->SetBinContent(fit+1,f->GetParameter(0));
+ fitNorm->SetBinError(fit+1,f->GetParError(0));
+ fitMu->SetBinContent(fit+1,f->GetParameter(1));
+ fitMu->SetBinError(fit+1,f->GetParError(1));
+ fitSigma->SetBinContent(fit+1,f->GetParameter(2));
+ fitSigma->SetBinError(fit+1,f->GetParError(2));
+ cout<<"\""<<f->GetName()<<"\","<<r[2]/r[3]<<","<<r[3]
+ <<","<<TMath::Prob(r[2],r[3])
+ <<","<<chi2/ndf
+ <<","<<ndf
+ <<","<<TMath::Prob(r[0],r[1])
+ <<","<<r[5];
+ for(int i=1;i<3;i++) {
+ cout<<","<<f->GetParameter(i)<<",\"\302\261\","<<f->GetParError(i);
+ }
+ cout<<"\n";
+ }
+
+ //=========================== L-curve ==========================
+ c4->cd(1);
+ lCurve->SetTitle("L curve;log_{10} L_{x};log_{10} L_{y}");
+ lCurve->SetLineColor(kRed);
+ lCurve->Draw("AL");
+ c4->cd(2);
+ gPad->Clear();
+ c4->cd(3);
+ logTauX->SetTitle(";log_{10} #tau;log_{10} L_{x}");
+ logTauX->SetLineColor(kBlue);
+ logTauX->Draw();
+ c4->cd(4);
+ logTauY->SetTitle(";log_{10} #tau;log_{10} L_{y}");
+ logTauY->SetLineColor(kBlue);
+ logTauY->Draw();
+ c4->SaveAs("lcurveL.eps");
+
+ //========================= rho and L-curve scan ===============
+ c4->cd(1);
+ logTauCurvature->SetTitle(";log_{10}(#tau);L curve curvature");
+ logTauCurvature->SetLineColor(kRed);
+ logTauCurvature->Draw();
+ c4->cd(2);
+ rhoScan->SetTitle(";log_{10}(#tau);average(#rho_{i})");
+ rhoScan->SetLineColor(kRed);
+ rhoScan->Draw();
+ c4->cd(3);
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
+ &table[iFitFALCurve].second);
+ c4->cd(4);
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFArhoO,"Tikhonov",tauFArho,
+ &table[iFitFArho].second);
+
+ c4->SaveAs("scanTau.eps");
+
+
+ TGraph *graphNiterAgoBay[4];
+ GetNiterGraphs(iAgoFirstFit,iAgoFirstFit+1,table,kRed-2,graphNiterAgoBay,20);
+ TGraph *graphNiterAgo[4];
+ GetNiterGraphs(iAgoFirstFit,iAgoLastFit,table,kRed,graphNiterAgo,24);
+#ifdef WITH_IDS
+ TGraph *graphNiterIDS[4];
+ GetNiterGraphs(iIDSFirstFit,iIDSLastFit,table,kMagenta,graphNiterIDS,21);
+#endif
+ TH1 *histNiterInversion[4];
+ GetNiterHist(iFitInversion,table,histNiterInversion,kCyan,1,1001);
+ TH1 *histNiterFALCurve[4];
+ GetNiterHist(iFitFALCurve,table,histNiterFALCurve,kBlue,2,3353);
+ TH1 *histNiterFArho[4];
+ GetNiterHist(iFitFArho,table,histNiterFArho,kAzure-4,3,3353);
+ TH1 *histNiterBinByBin[4];
+ GetNiterHist(iFitBinByBin,table,histNiterBinByBin,kOrange+1,3,3335);
+
+ histNiterInversion[0]->GetYaxis()->SetRangeUser(0.3,500.);
+ histNiterInversion[1]->GetYaxis()->SetRangeUser(-0.1,0.9);
+ histNiterInversion[2]->GetYaxis()->SetRangeUser(5.6,6.3);
+ histNiterInversion[3]->GetYaxis()->SetRangeUser(1.6,2.4);
+
+ TLine *line=0;
+ c1->cd();
+ for(int i=0;i<2;i++) {
+ gPad->Clear();
+ gPad->SetLogx();
+ gPad->SetLogy((i<1));
+ if(! histNiterInversion[i]) continue;
+ histNiterInversion[i]->Draw("][");
+ histNiterFALCurve[i]->Draw("SAME ][");
+ histNiterFArho[i]->Draw("SAME ][");
+ histNiterBinByBin[i]->Draw("SAME ][");
+ graphNiterAgo[i]->Draw("LP");
+ graphNiterAgoBay[i]->SetMarkerStyle(20);
+ graphNiterAgoBay[i]->Draw("P");
+#ifdef WITH_IDS
+ graphNiterIDS[i]->Draw("LP");
+#endif
+ TLegend *legend;
+ if(i==1) {
+ legend=new TLegend(0.48,0.28,0.87,0.63);
+ } else {
+ legend=new TLegend(0.45,0.5,0.88,0.88);
+ }
+ legend->SetBorderSize(0);
+ legend->SetFillStyle(1001);
+ legend->SetFillColor(kWhite);
+ legend->SetTextSize(kLegendFontSize*0.7);
+ legend->AddEntry( histNiterInversion[0],"inversion","l");
+ legend->AddEntry( histNiterFALCurve[0],"Tikhonov L-curve","l");
+ legend->AddEntry( histNiterFArho[0],"Tikhonov global cor.","l");
+ legend->AddEntry( histNiterBinByBin[0],"bin-by-bin","l");
+ legend->AddEntry( graphNiterAgoBay[0],"\"Bayesian\"","p");
+ legend->AddEntry( graphNiterAgo[0],"iterative","p");
+#ifdef WITH_IDS
+ legend->AddEntry( graphNiterIDS[0],"IDS","p");
+#endif
+ legend->Draw();
+
+ c1->SaveAs(TString::Format("niter%d.eps",i));
+ }
+
+ //c4->cd(1);
+ //DrawPadCorrelations(histRhoFALCurveO,&table);
+ c2sq->cd(1);
+ DrawPadTruth(histMcsigGenO,histDataGenO,histOutputFALCurveO,"Tikhonov",tauFA,
+ &table[iFitFALCurve].second,table[iFitFALCurve].first);
+
+ c2sq->cd(2);
+ gPad->SetLogx();
+ gPad->SetLogy(false);
+ histNiterInversion[3]->DrawClone("E2");
+ histNiterInversion[3]->SetFillStyle(0);
+ histNiterInversion[3]->Draw("SAME HIST ][");
+ histNiterFALCurve[3]->DrawClone("SAME E2");
+ histNiterFALCurve[3]->SetFillStyle(0);
+ histNiterFALCurve[3]->Draw("SAME HIST ][");
+ histNiterFArho[3]->DrawClone("SAME E2");
+ histNiterFArho[3]->SetFillStyle(0);
+ histNiterFArho[3]->Draw("SAME HIST ][");
+ histNiterBinByBin[3]->DrawClone("SAME E2");
+ histNiterBinByBin[3]->SetFillStyle(0);
+ histNiterBinByBin[3]->Draw("SAME HIST ][");
+ double yTrue=1.8;
+ line=new TLine(histNiterInversion[3]->GetXaxis()->GetXmin(),
+ yTrue,
+ histNiterInversion[3]->GetXaxis()->GetXmax(),
+ yTrue);
+ line->SetLineWidth(3);
+ line->Draw();
+
+ graphNiterAgo[3]->Draw("LP");
+ graphNiterAgoBay[3]->SetMarkerStyle(20);
+ graphNiterAgoBay[3]->Draw("P");
+#ifdef WITH_IDS
+ graphNiterIDS[3]->Draw("LP");
+#endif
+
+ TLegend *legend;
+ legend=new TLegend(0.55,0.53,0.95,0.97);
+ legend->SetBorderSize(0);
+ legend->SetFillStyle(1001);
+ legend->SetFillColor(kWhite);
+ legend->SetTextSize(kLegendFontSize);
+ legend->AddEntry( line,"truth","l");
+ legend->AddEntry( histNiterInversion[3],"inversion","l");
+ legend->AddEntry( histNiterFALCurve[3],"Tikhonov L-curve","l");
+ legend->AddEntry( histNiterFArho[3],"Tikhonov global cor.","l");
+ legend->AddEntry( histNiterBinByBin[3],"bin-by-bin","l");
+ legend->AddEntry( graphNiterAgoBay[3],"\"Bayesian\"","p");
+ legend->AddEntry( graphNiterAgo[3],"iterative","p");
+#ifdef WITH_IDS
+ legend->AddEntry( graphNiterIDS[3],"IDS","p");
+#endif
+ legend->Draw();
+
+ c2sq->SaveAs("fitSigma.eps");
+
+ outputFile->Write();
+
+ delete outputFile;
+
+}
+
+void GetNiterGraphs(int iFirst,int iLast,vector<pair<TF1*,
+ vector<double> > > const &table,int color,
+ TGraph *graph[4],int style) {
+ TVectorD niter(iLast-iFirst);
+ TVectorD eniter(iLast-iFirst);
+ TVectorD chi2(iLast-iFirst);
+ TVectorD gcor(iLast-iFirst);
+ TVectorD mean(iLast-iFirst);
+ TVectorD emean(iLast-iFirst);
+ TVectorD sigma(iLast-iFirst);
+ TVectorD esigma(iLast-iFirst);
+ for(int ifit=iFirst;ifit<iLast;ifit++) {
+ vector<double> const &r=table[ifit].second;
+ niter(ifit-iFirst)=r[4];
+ chi2(ifit-iFirst)=r[2]/r[3];
+ gcor(ifit-iFirst)=r[5];
+ TF1 const *f=table[ifit].first;
+ mean(ifit-iFirst)=f->GetParameter(1);
+ emean(ifit-iFirst)=f->GetParError(1);
+ sigma(ifit-iFirst)=f->GetParameter(2);
+ esigma(ifit-iFirst)=f->GetParError(2);
+ }
+ graph[0]=new TGraph(niter,chi2);
+ graph[1]=new TGraph(niter,gcor);
+ graph[2]=new TGraphErrors(niter,mean,eniter,emean);
+ graph[3]=new TGraphErrors(niter,sigma,eniter,esigma);
+ for(int g=0;g<4;g++) {
+ if(graph[g]) {
+ graph[g]->SetLineColor(color);
+ graph[g]->SetMarkerColor(color);
+ graph[g]->SetMarkerStyle(style);
+ }
+ }
+}
+
+void GetNiterHist(int ifit,vector<pair<TF1*,vector<double> > > const &table,
+ TH1 *hist[4],int color,int style,int fillStyle) {
+ vector<double> const &r=table[ifit].second;
+ TF1 const *f=table[ifit].first;
+ hist[0]=new TH1D(table[ifit].first->GetName()+TString("_chi2"),
+ ";iteration;unfold-truth #chi^{2}/N_{D.F.}",1,0.2,1500.);
+ hist[0]->SetBinContent(1,r[2]/r[3]);
+
+ hist[1]=new TH1D(table[ifit].first->GetName()+TString("_gcor"),
+ ";iteration;avg(#rho_{i})",1,0.2,1500.);
+ hist[1]->SetBinContent(1,r[5]);
+ hist[2]=new TH1D(table[ifit].first->GetName()+TString("_mu"),
+ ";iteration;parameter #mu",1,0.2,1500.);
+ hist[2]->SetBinContent(1,f->GetParameter(1));
+ hist[2]->SetBinError(1,f->GetParError(1));
+ hist[3]=new TH1D(table[ifit].first->GetName()+TString("_sigma"),
+ ";iteration;parameter #sigma",1,0.2,1500.);
+ hist[3]->SetBinContent(1,f->GetParameter(2));
+ hist[3]->SetBinError(1,f->GetParError(2));
+ for(int h=0;h<4;h++) {
+ if(hist[h]) {
+ hist[h]->SetLineColor(color);
+ hist[h]->SetLineStyle(style);
+ if( hist[h]->GetBinError(1)>0.0) {
+ hist[h]->SetFillColor(color-10);
+ hist[h]->SetFillStyle(fillStyle);
+ }
+ hist[h]->SetMarkerStyle(0);
+ }
+ }
+}
+
+void CreateHistogramCopies(TH1 *h[3],TUnfoldBinning const *binning) {
+ TString baseName(h[0]->GetName());
+ Int_t *binMap;
+ h[1]=binning->CreateHistogram(baseName+"_axis",kTRUE,&binMap);
+ h[2]=(TH1 *)h[1]->Clone(baseName+"_binw");
+ Int_t nMax=binning->GetEndBin()+1;
+ for(Int_t iSrc=0;iSrc<nMax;iSrc++) {
+ Int_t iDest=binMap[iSrc];
+ double c=h[0]->GetBinContent(iSrc)+h[1]->GetBinContent(iDest);
+ double e=TMath::Hypot(h[0]->GetBinError(iSrc),h[1]->GetBinError(iDest));
+ h[1]->SetBinContent(iDest,c);
+ h[1]->SetBinError(iDest,e);
+ h[2]->SetBinContent(iDest,c);
+ h[2]->SetBinError(iDest,e);
+ }
+ for(int iDest=0;iDest<=h[2]->GetNbinsX()+1;iDest++) {
+ double c=h[2]->GetBinContent(iDest);
+ double e=h[2]->GetBinError(iDest);
+ double bw=binning->GetBinSize(iDest);
+ /* if(bw!=h[2]->GetBinWidth(iDest)) {
+ cout<<"bin "<<iDest<<" width="<<bw<<" "<<h[2]->GetBinWidth(iDest)
+ <<"\n";
+ } */
+ if(bw>0.0) {
+ h[2]->SetBinContent(iDest,c/bw);
+ h[2]->SetBinError(iDest,e/bw);
+ } else {
+ }
+ }
+}
+
+void CreateHistogramCopies(TH2 *h[3],TUnfoldBinning const *binningX) {
+ if(binningX) {
+ h[2]=0;
+ h[3]=0;
+ } else {
+ h[2]=0;
+ h[3]=0;
+ }
+}
+
+TH2 *AddOverflowXY(TH2 *h,double widthX,double widthY) {
+ // add overflow bin to X-axis
+ int nx=h->GetNbinsX();
+ int ny=h->GetNbinsY();
+ double *xBins=new double[nx+2];
+ double *yBins=new double[ny+2];
+ for(int i=1;i<=nx;i++) {
+ xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
+ }
+ xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
+ xBins[nx+1]=xBins[nx]+widthX;
+ for(int i=1;i<=ny;i++) {
+ yBins[i-1]=h->GetYaxis()->GetBinLowEdge(i);
+ }
+ yBins[ny]=h->GetYaxis()->GetBinUpEdge(ny);
+ yBins[ny+1]=yBins[ny]+widthY;
+ TString name(h->GetName());
+ name+="U";
+ TH2 *r=new TH2D(name,h->GetTitle(),nx+1,xBins,ny+1,yBins);
+ for(int ix=0;ix<=nx+1;ix++) {
+ for(int iy=0;iy<=ny+1;iy++) {
+ r->SetBinContent(ix,iy,h->GetBinContent(ix,iy));
+ r->SetBinError(ix,iy,h->GetBinError(ix,iy));
+ }
+ }
+ delete [] yBins;
+ delete [] xBins;
+ return r;
+}
+
+TH1 *AddOverflowX(TH1 *h,double widthX) {
+ // add overflow bin to X-axis
+ int nx=h->GetNbinsX();
+ double *xBins=new double[nx+2];
+ for(int i=1;i<=nx;i++) {
+ xBins[i-1]=h->GetXaxis()->GetBinLowEdge(i);
+ }
+ xBins[nx]=h->GetXaxis()->GetBinUpEdge(nx);
+ xBins[nx+1]=xBins[nx]+widthX;
+ TString name(h->GetName());
+ name+="U";
+ TH1 *r=new TH1D(name,h->GetTitle(),nx+1,xBins);
+ for(int ix=0;ix<=nx+1;ix++) {
+ r->SetBinContent(ix,h->GetBinContent(ix));
+ r->SetBinError(ix,h->GetBinError(ix));
+ }
+ delete [] xBins;
+ return r;
+}
+
+void DrawOverflowX(TH1 *h,double posy) {
+ double x1=h->GetXaxis()->GetBinLowEdge(h->GetNbinsX());
+ double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
+ double y0=h->GetYaxis()->GetBinLowEdge(1);
+ double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
+ if(h->GetDimension()==1) {
+ y0=h->GetMinimum();
+ y2=h->GetMaximum();
+ }
+ double w1=-0.3;
+ TText *textX=new TText((1.+w1)*x2-w1*x1,(1.-posy)*y0+posy*y2,"Overflow bin");
+ textX->SetNDC(kFALSE);
+ textX->SetTextSize(0.05);
+ textX->SetTextAngle(90.);
+ textX->Draw();
+ TLine *lineX=new TLine(x1,y0,x1,y2);
+ lineX->Draw();
+}
+
+void DrawOverflowY(TH1 *h,double posx) {
+ double x0=h->GetXaxis()->GetBinLowEdge(1);
+ double x2=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
+ double y1=h->GetYaxis()->GetBinLowEdge(h->GetNbinsY());;
+ double y2=h->GetYaxis()->GetBinUpEdge(h->GetNbinsY());;
+ double w1=-0.3;
+ TText *textY=new TText((1.-posx)*x0+posx*x2,(1.+w1)*y1-w1*y2,"Overflow bin");
+ textY->SetNDC(kFALSE);
+ textY->SetTextSize(0.05);
+ textY->Draw();
+ TLine *lineY=new TLine(x0,y1,x2,y1);
+ lineY->Draw();
+}
+
+void DrawPadProbability(TH2 *h) {
+ h->Draw("COLZ");
+ h->SetTitle("migration probabilities;P_{T}(gen) [GeV];P_{T}(rec) [GeV]");
+ DrawOverflowX(h,0.05);
+ DrawOverflowY(h,0.35);
+}
+
+void DrawPadEfficiency(TH1 *h) {
+ h->SetTitle("efficiency;P_{T}(gen) [GeV];#epsilon");
+ h->SetLineColor(kBlue);
+ h->SetMinimum(0.75);
+ h->SetMaximum(1.0);
+ h->Draw();
+ DrawOverflowX(h,0.05);
+ TLegend *legEfficiency=new TLegend(0.3,0.58,0.6,0.75);
+ legEfficiency->SetBorderSize(0);
+ legEfficiency->SetFillStyle(0);
+ legEfficiency->SetTextSize(kLegendFontSize);
+ legEfficiency->AddEntry(h,"reconstruction","l");
+ legEfficiency->AddEntry((TObject*)0," efficiency","");
+ legEfficiency->Draw();
+}
+
+void DrawPadReco(TH1 *histMcRec,TH1 *histMcbgrRec,TH1 *histDataRec,
+ TH1 *histDataUnfold,TH2 *histProbability,TH2 *histRhoij) {
+ //gPad->SetLogy(kTRUE);
+ double amax=0.0;
+ for(int i=1;i<=histMcRec->GetNbinsX();i++) {
+ amax=TMath::Max(amax,histMcRec->GetBinContent(i)
+ +5.0*histMcRec->GetBinError(i));
+ amax=TMath::Max(amax,histDataRec->GetBinContent(i)
+ +2.0*histDataRec->GetBinError(i));
+ }
+ histMcRec->SetTitle("Reconstructed;P_{T}(rec);Nevent / GeV");
+ histMcRec->Draw("HIST");
+ histMcRec->SetLineColor(kBlue);
+ histMcRec->SetMinimum(1.0);
+ histMcRec->SetMaximum(amax);
+ //histMcbgrRec->SetFillMode(1);
+ histMcbgrRec->SetLineColor(kBlue-6);
+ histMcbgrRec->SetFillColor(kBlue-10);
+ histMcbgrRec->Draw("SAME HIST");
+
+ TH1 * histFoldBack=0;
+ if(histDataUnfold && histProbability && histRhoij) {
+ histFoldBack=(TH1 *)
+ histMcRec->Clone(histDataUnfold->GetName()+TString("_folded"));
+ int nrec=histFoldBack->GetNbinsX();
+ if((nrec==histProbability->GetNbinsY())&&
+ (nrec==histMcbgrRec->GetNbinsX())&&
+ (nrec==histDataRec->GetNbinsX())
+ ) {
+ for(int ix=1;ix<=nrec;ix++) {
+ double sum=0.0;
+ double sume2=0.0;
+ for(int iy=0;iy<=histProbability->GetNbinsX()+1;iy++) {
+ sum += histDataUnfold->GetBinContent(iy)*
+ histDataUnfold->GetBinWidth(iy)*
+ histProbability->GetBinContent(iy,ix);
+ for(int iy2=0;iy2<=histProbability->GetNbinsX()+1;iy2++) {
+ sume2 += histDataUnfold->GetBinError(iy)*
+ histDataUnfold->GetBinWidth(iy)*
+ histProbability->GetBinContent(iy,ix)*
+ histDataUnfold->GetBinError(iy2)*
+ histDataUnfold->GetBinWidth(iy2)*
+ histProbability->GetBinContent(iy2,ix)*
+ histRhoij->GetBinContent(iy,iy2);
+ }
+ }
+ sum /= histFoldBack->GetBinWidth(ix);
+ sum += histMcbgrRec->GetBinContent(ix);
+ histFoldBack->SetBinContent(ix,sum);
+ histFoldBack->SetBinError(ix,TMath::Sqrt(sume2)
+ /histFoldBack->GetBinWidth(ix));
+ }
+ } else {
+ cout<<"can not fold back: "<<nrec
+ <<" "<<histProbability->GetNbinsY()
+ <<" "<<histMcbgrRec->GetNbinsX()
+ <<" "<<histDataRec->GetNbinsX()
+ <<"\n";
+ exit(0);
+ }
+
+ histFoldBack->SetLineColor(kBlack);
+ histFoldBack->SetMarkerStyle(0);
+ histFoldBack->Draw("SAME HIST");
+ }
+
+ histDataRec->SetLineColor(kRed);
+ histDataRec->SetMarkerColor(kRed);
+ histDataRec->Draw("SAME");
+ DrawOverflowX(histMcRec,0.5);
+
+ TLegend *legRec=new TLegend(0.4,0.5,0.68,0.85);
+ legRec->SetBorderSize(0);
+ legRec->SetFillStyle(0);
+ legRec->SetTextSize(kLegendFontSize);
+ legRec->AddEntry(histMcRec,"MC total","l");
+ legRec->AddEntry(histMcbgrRec,"background","f");
+ if(histFoldBack) {
+ int ndf=-kNbinC;
+ double sumD=0.,sumF=0.,chi2=0.;
+ for(int i=1;i<=histDataRec->GetNbinsX();i++) {
+ //cout<<histDataRec->GetBinContent(i)<<" "<<histFoldBack->GetBinContent(i)<<" "<<" w="<<histFoldBack->GetBinWidth(i)<<"\n";
+ sumD+=histDataRec->GetBinContent(i)*histDataRec->GetBinWidth(i);
+ sumF+=histFoldBack->GetBinContent(i)*histFoldBack->GetBinWidth(i);
+ double pull=(histFoldBack->GetBinContent(i)-histDataRec->GetBinContent(i))/histDataRec->GetBinError(i);
+ chi2+= pull*pull;
+ ndf+=1;
+ }
+ legRec->AddEntry(histDataRec,TString::Format("data N_{evt}=%.0f",sumD),"lp");
+ legRec->AddEntry(histFoldBack,TString::Format("folded N_{evt}=%.0f",sumF),"l");
+ legRec->AddEntry((TObject*)0,TString::Format("#chi^{2}=%.1f ndf=%d",chi2,ndf),"");
+ //exit(0);
+ } else {
+ legRec->AddEntry(histDataRec,"data","lp");
+ }
+ legRec->Draw();
+}
+
+void DrawPadTruth(TH1 *histMcsigGen,TH1 *histDataGen,TH1 *histDataUnfold,
+ char const *text,double tau,vector<double> const *r,
+ TF1 *f) {
+ //gPad->SetLogy(kTRUE);
+ double amin=0.;
+ double amax=0.;
+ for(int i=1;i<=histMcsigGen->GetNbinsX();i++) {
+ if(histDataUnfold) {
+ amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
+ -1.1*histDataUnfold->GetBinError(i));
+ amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
+ +1.1*histDataUnfold->GetBinError(i));
+ }
+ amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
+ -histMcsigGen->GetBinError(i));
+ amin=TMath::Min(amin,histDataGen->GetBinContent(i)
+ -histDataGen->GetBinError(i));
+ amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
+ +10.*histMcsigGen->GetBinError(i));
+ amax=TMath::Max(amax,histDataGen->GetBinContent(i)
+ +2.*histDataGen->GetBinError(i));
+ }
+ histMcsigGen->SetMinimum(amin);
+ histMcsigGen->SetMaximum(amax);
+
+ histMcsigGen->SetTitle("Truth;P_{T};Nevent / GeV");
+ histMcsigGen->SetLineColor(kBlue);
+ histMcsigGen->Draw("HIST");
+ histDataGen->SetLineColor(kRed);
+ histDataGen->SetMarkerColor(kRed);
+ histDataGen->SetMarkerSize(1.0);
+ histDataGen->Draw("SAME HIST");
+ if(histDataUnfold) {
+ histDataUnfold->SetMarkerStyle(21);
+ histDataUnfold->SetMarkerSize(0.7);
+ histDataUnfold->Draw("SAME");
+ }
+ DrawOverflowX(histMcsigGen,0.5);
+
+ if(f) {
+ f->SetLineStyle(1);
+ f->Draw("SAME");
+ }
+
+ TLegend *legTruth=new TLegend(0.32,0.65,0.6,0.9);
+ legTruth->SetBorderSize(0);
+ legTruth->SetFillStyle(0);
+ legTruth->SetTextSize(kLegendFontSize);
+ legTruth->AddEntry(histMcsigGen,"MC","l");
+ if(!histDataUnfold) legTruth->AddEntry((TObject *)0," Landau(5,2)","");
+ legTruth->AddEntry(histDataGen,"data","l");
+ if(!histDataUnfold) legTruth->AddEntry((TObject *)0," Landau(6,1.8)","");
+ if(histDataUnfold) {
+ TString t;
+ if(text) t=text;
+ else t=histDataUnfold->GetName();
+ if(tau>0) {
+ t+=TString::Format(" #tau=%.2g",tau);
+ }
+ legTruth->AddEntry(histDataUnfold,t,"lp");
+ if(r) {
+ legTruth->AddEntry((TObject *)0,"test wrt data:","");
+ legTruth->AddEntry((TObject *)0,TString::Format
+ ("#chi^{2}/%d=%.1f prob=%.3f",
+ (int)(*r)[3],(*r)[2]/(*r)[3],
+ TMath::Prob((*r)[2],(*r)[3])),"");
+ }
+ }
+ if(f) {
+ legTruth->AddEntry(f,"fit","l");
+ }
+ legTruth->Draw();
+ if(histDataUnfold ) {
+ TPad *subpad = new TPad("subpad","",0.35,0.29,0.88,0.68);
+ subpad->SetFillStyle(0);
+ subpad->Draw();
+ subpad->cd();
+ amin=0.;
+ amax=0.;
+ int istart=11;
+ for(int i=istart;i<=histMcsigGen->GetNbinsX();i++) {
+ amin=TMath::Min(amin,histMcsigGen->GetBinContent(i)
+ -histMcsigGen->GetBinError(i));
+ amin=TMath::Min(amin,histDataGen->GetBinContent(i)
+ -histDataGen->GetBinError(i));
+ amin=TMath::Min(amin,histDataUnfold->GetBinContent(i)
+ -histDataUnfold->GetBinError(i));
+ amax=TMath::Max(amax,histMcsigGen->GetBinContent(i)
+ +histMcsigGen->GetBinError(i));
+ amax=TMath::Max(amax,histDataGen->GetBinContent(i)
+ +histDataGen->GetBinError(i));
+ amax=TMath::Max(amax,histDataUnfold->GetBinContent(i)
+ +histDataUnfold->GetBinError(i));
+ }
+ TH1 *copyMcsigGen=(TH1*)histMcsigGen->Clone();
+ TH1 *copyDataGen=(TH1*)histDataGen->Clone();
+ TH1 *copyDataUnfold=(TH1*)histDataUnfold->Clone();
+ copyMcsigGen->GetXaxis()->SetRangeUser
+ (copyMcsigGen->GetXaxis()->GetBinLowEdge(istart),
+ copyMcsigGen->GetXaxis()->GetBinUpEdge(copyMcsigGen->GetNbinsX()-1));
+ copyMcsigGen->SetTitle(";;");
+ copyMcsigGen->GetYaxis()->SetRangeUser(amin,amax);
+ copyMcsigGen->Draw("HIST");
+ copyDataGen->Draw("SAME HIST");
+ copyDataUnfold->Draw("SAME");
+ if(f) {
+ ((TF1 *)f->Clone())->Draw("SAME");
+ }
+ }
+}
+
+void DrawPadCorrelations(TH2 *h,
+ vector<pair<TF1*,vector<double> > > const *table) {
+ h->SetMinimum(-1.);
+ h->SetMaximum(1.);
+ h->SetTitle("correlation coefficients;P_{T}(gen) [GeV];P_{T}(gen) [GeV]");
+ h->Draw("COLZ");
+ DrawOverflowX(h,0.05);
+ DrawOverflowY(h,0.05);
+ if(table) {
+ TLegend *legGCor=new TLegend(0.13,0.6,0.5,0.8);
+ legGCor->SetBorderSize(0);
+ legGCor->SetFillStyle(0);
+ legGCor->SetTextSize(kLegendFontSize);
+ vector<double> const &r=(*table)[table->size()-1].second;
+ legGCor->AddEntry((TObject *)0,TString::Format("min(#rho_{ij})=%5.2f",r[6]),"");
+ legGCor->AddEntry((TObject *)0,TString::Format("max(#rho_{ij})=%5.2f",r[7]),"");
+ legGCor->AddEntry((TObject *)0,TString::Format("avg(#rho_i)=%5.2f",r[5]),"");
+ legGCor->Draw();
+ }
+}
+
+TH1 *g_fcnHist=0;
+TMatrixD *g_fcnMatrix=0;
+
+void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag) {
+ if(flag==0) {
+ cout<<"fcn flag=0: npar="<<npar<<" gin="<<gin<<" par=[";
+ for(int i=0;i<npar;i++) {
+ cout<<" "<<u[i];
+ }
+ cout<<"]\n";
+ }
+ int n=g_fcnMatrix->GetNrows();
+ TVectorD dy(n);
+ double x0=0,y0=0.;
+ for(int i=0;i<=n;i++) {
+ double x1;
+ if(i<1) x1=g_fcnHist->GetXaxis()->GetBinLowEdge(i+1);
+ else x1=g_fcnHist->GetXaxis()->GetBinUpEdge(i);
+ double y1=TMath::LandauI((x1-u[1])/u[2]);
+ if(i>0) {
+ double iy=u[0]*u[2]*(y1-y0)/(x1-x0);
+ dy(i-1)=iy-g_fcnHist->GetBinContent(i);
+ //cout<<"i="<<i<<" iy="<<iy<<" delta="<< dy(i-1)<<"\n";
+ }
+ x0=x1;
+ y0=y1;
+ //cout<<"i="<<i<<" y1="<<y1<<" x1="<<x1<<"\n";
+ }
+ TVectorD Hdy=(*g_fcnMatrix) * dy;
+ //Hdy.Print();
+ f=Hdy*dy;
+ //exit(0);
+}
+
+
+TFitResultPtr DoFit(TH1 *h,TH2 *rho,TH1 *truth,const char *text,
+ vector<pair<TF1 *,vector<double> > > &table,int niter) {
+ TString option="IESN";
+ cout<<h->GetName()<<"\n";
+ double gcorAvg=0.;
+ double rhoMin=0.;
+ double rhoMax=0.;
+ if(rho) {
+ g_fcnHist=h;
+ int n=h->GetNbinsX()-1; // overflow is included as extra bin, exclude in fit
+ TMatrixDSym v(n);
+ //g_fcnMatrix=new TMatrixD(n,n);
+ for(int i=0;i<n;i++) {
+ for(int j=0;j<n;j++) {
+ v(i,j)=rho->GetBinContent(i+1,j+1)*
+ (h->GetBinError(i+1)*h->GetBinError(j+1));
+ }
+ }
+ TMatrixDSymEigen ev(v);
+ TMatrixD d(n,n);
+ TVectorD di(ev.GetEigenValues());
+ for(int i=0;i<n;i++) {
+ if(di(i)>0.0) {
+ d(i,i)=1./di(i);
+ } else {
+ cout<<"bad eigenvalue i="<<i<<" di="<<di(i)<<"\n";
+ exit(0);
+ }
+ }
+ TMatrixD O(ev.GetEigenVectors());
+ TMatrixD DOT(d,TMatrixD::kMultTranspose,O);
+ g_fcnMatrix=new TMatrixD(O,TMatrixD::kMult,DOT);
+ TMatrixD test(*g_fcnMatrix,TMatrixD::kMult,v);
+ int error=0;
+ for(int i=0;i<n;i++) {
+ if(TMath::Abs(test(i,i)-1.0)>1.E-7) {
+ error++;
+ }
+ for(int j=0;j<n;j++) {
+ if(i==j) continue;
+ if(TMath::Abs(test(i,j)>1.E-7)) error++;
+ }
+ }
+ // calculate global correlation coefficient (all bins)
+ TMatrixDSym v1(n+1);
+ rhoMin=1.;
+ rhoMax=-1.;
+ for(int i=0;i<=n;i++) {
+ for(int j=0;j<=n;j++) {
+ double rho_ij=rho->GetBinContent(i+1,j+1);
+ v1(i,j)=rho_ij*
+ (h->GetBinError(i+1)*h->GetBinError(j+1));
+ if(i!=j) {
+ if(rho_ij<rhoMin) rhoMin=rho_ij;
+ if(rho_ij>rhoMax) rhoMax=rho_ij;
+ }
+ }
+ }
+ TMatrixDSymEigen ev1(v1);
+ TMatrixD d1(n+1,n+1);
+ TVectorD di1(ev1.GetEigenValues());
+ for(int i=0;i<=n;i++) {
+ if(di1(i)>0.0) {
+ d1(i,i)=1./di1(i);
+ } else {
+ cout<<"bad eigenvalue i="<<i<<" di1="<<di1(i)<<"\n";
+ exit(0);
+ }
+ }
+ TMatrixD O1(ev1.GetEigenVectors());
+ TMatrixD DOT1(d1,TMatrixD::kMultTranspose,O1);
+ TMatrixD vinv1(O1,TMatrixD::kMult,DOT1);
+ for(int i=0;i<=n;i++) {
+ double gcor2=1.-1./(vinv1(i,i)*v1(i,i));
+ if(gcor2>=0.0) {
+ double gcor=TMath::Sqrt(gcor2);
+ gcorAvg += gcor;
+ } else {
+ cout<<"bad global correlation "<<i<<" "<<gcor2<<"\n";
+ }
+ }
+ gcorAvg /=(n+1);
+ /* if(error) {
+ v.Print();
+ g_fcnMatrix->Print();
+ exit(0);
+ } */
+ //g_fcnMatrix->Invert();
+ //from: HFitImpl.cxx
+ // TVirtualFitter::FCNFunc_t userFcn = 0;
+ // typedef void (* FCNFunc_t )(Int_t &npar, Double_t *gin, Double_t &f, Double_t *u, Int_t flag);
+ // userFcn = (TVirtualFitter::GetFitter())->GetFCN();
+ // (TVirtualFitter::GetFitter())->SetUserFunc(f1);
+ //...
+ //fitok = fitter->FitFCN( userFcn );
+
+ TVirtualFitter::Fitter(h)->SetFCN(fcn);
+ option += "U";
+
+ }
+ double xmax=h->GetXaxis()->GetBinUpEdge(h->GetNbinsX()-1);
+ TF1 *landau=new TF1(text,"[0]*TMath::Landau(x,[1],[2],0)",
+ 0.,xmax);
+ landau->SetParameter(0,6000.);
+ landau->SetParameter(1,5.);
+ landau->SetParameter(2,2.);
+ landau->SetParError(0,10.);
+ landau->SetParError(1,0.5);
+ landau->SetParError(2,0.1);
+ TFitResultPtr s=h->Fit(landau,option,0,0.,xmax);
+ vector<double> r(8);
+ int np=landau->GetNpar();
+ fcn(np,0,r[0],landau->GetParameters(),0);
+ r[1]=h->GetNbinsX()-1-landau->GetNpar();
+ for(int i=0;i<h->GetNbinsX()-1;i++) {
+ double di=h->GetBinContent(i+1)-truth->GetBinContent(i+1);
+ if(g_fcnMatrix) {
+ for(int j=0;j<h->GetNbinsX()-1;j++) {
+ double dj=h->GetBinContent(j+1)-truth->GetBinContent(j+1);
+ r[2]+=di*dj*(*g_fcnMatrix)(i,j);
+ }
+ } else {
+ double pull=di/h->GetBinError(i+1);
+ r[2]+=pull*pull;
+ }
+ r[3]+=1.0;
+ }
+ r[4]=niter;
+ if(!niter) r[4]=0.25;
+ r[5]=gcorAvg;
+ r[6]=rhoMin;
+ r[7]=rhoMax;
+ if(rho) {
+ g_fcnHist=0;
+ delete g_fcnMatrix;
+ g_fcnMatrix=0;
+ }
+ table.push_back(make_pair(landau,r));
+ return s;
+}
+
+#ifdef WITH_IDS
+
+//===================== interface to IDS unfolding code follows here
+// contact Bogdan Malescu to find it
+
+#include "ids_code.cc"
+
+void IDSfirst(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, Double_t lambdaL_, TVectorD* &unfres1IDS_,TVectorD *&soustr){
+
+ int N_=data->GetNrows();
+ soustr = new TVectorD(N_);
+ for( Int_t i=0; i<N_; i++ ){ (*soustr)[i] = 0.; }
+ unfres1IDS_ = Unfold( data, dataErr, A_, N_, lambdaL_, soustr );
+}
+
+void IDSiterate(TVectorD *data, TVectorD *dataErr, TMatrixD *A_, TMatrixD *Am_, Double_t lambdaU_, Double_t lambdaM_, Double_t lambdaS_,TVectorD* &unfres2IDS_ ,TVectorD *&soustr) {
+
+ int N_=data->GetNrows();
+ ModifyMatrix( Am_, A_, unfres2IDS_, dataErr, N_, lambdaM_, soustr, lambdaS_ );
+ delete unfres2IDS_;
+ unfres2IDS_ = Unfold( data, dataErr, Am_, N_, lambdaU_, soustr );
+}
+
+#endif
--
GitLab