diff --git a/tutorials/math/testUnfold1.C b/tutorials/math/testUnfold1.C
deleted file mode 100644
index a76dec9ea632d19ef6d429c413fe7e85e03c800a..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold1.C
+++ /dev/null
@@ -1,496 +0,0 @@
-/// \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 decsribing 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: unfoldede 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 log(tau)
-/// - the star indicates the final choice of tau
-/// 6. the L curve
-///
-/// Version 17.0, updated for using the classes TUnfoldDensity, TUnfoldBinning
-///
-/// History:
-/// - 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
-///
-/// \macro_image
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#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
- // Corelated 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 withd 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 scame
- // 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);
-
- //-------------------------------------------------------
- // retreive 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/math/testUnfold2.C b/tutorials/math/testUnfold2.C
deleted file mode 100644
index fcdb22751a72244a94bbf793a4bfcc7ceed839f6..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold2.C
+++ /dev/null
@@ -1,344 +0,0 @@
-/// \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
-///
-/// Version 17.0, updated for changed methods in TUnfold
-///
-/// History:
-/// - Version 16.1, parallel to changes in TUnfold
-/// - Version 16.0, parallel to changes in TUnfold
-/// - Version 15, with automatic L-curve scan, simplified example
-/// - Version 14, with changes in TUnfoldSys.cxx
-/// - Version 13, with changes to TUnfold.C
-/// - Version 12, with improvements to TUnfold.cxx
-/// - Version 11, print chi**2 and number of degrees of freedom
-/// - Version 10, with bug-fix in TUnfold.cxx
-/// - Version 9, with bug-fix in TUnfold.cxx, TUnfold.h
-/// - Version 8, with bug-fix in TUnfold.cxx, TUnfold.h
-/// - Version 7, with bug-fix in TUnfold.cxx, TUnfold.h
-/// - Version 6a, fix problem with dynamic array allocation under windows
-/// - Version 6, re-include class MyUnfold in the example
-/// - Version 5, move class MyUnfold to seperate files
-/// - Version 4, with bug-fix in TUnfold.C
-/// - Version 3, with bug-fix in TUnfold.C
-/// - Version 2, with changed ScanLcurve() arguments
-/// - Version 1, remove L curve analysis, use ScanLcurve() method instead
-/// - Version 0, L curve analysis included here
-///
-/// \macro_image
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#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 correspoinding to estimated peak position
- Int_t iPeek=(Int_t)(nGen*(estimatedPeakPosition-xminGen)/(xmaxGen-xminGen)
- // offset 1.5
- // accounts for start bin 1
- // and rounding errors +0.5
- +1.5);
- // regularize output bins 1..iPeek-nPeek
- unfold.RegularizeBins(1,1,iPeek-nPeek,regMode);
- // regularize output bins iPeek+nPeek..nGen
- unfold.RegularizeBins(iPeek+nPeek,1,nGen-(iPeek+nPeek),regMode);
-
- // unfolding
- //-----------
-
- // set input distribution and bias scale (=0)
- if(unfold.SetInput(histMdetData,0.0)>=10000) {
- std::cout<<"Unfolding result may be wrong\n";
- }
-
- // do the unfolding here
- Double_t tauMin=0.0;
- Double_t tauMax=0.0;
- Int_t nScan=30;
- Int_t iBest;
- TSpline *logTauX,*logTauY;
- TGraph *lCurve;
- // this method scans the parameter tau and finds the kink in the L curve
- // finally, the unfolding is done for the "best" choice of tau
- iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
- std::cout<<"tau="<<unfold.GetTau()<<"\n";
- std::cout<<"chi**2="<<unfold.GetChi2A()<<"+"<<unfold.GetChi2L()
- <<" / "<<unfold.GetNdf()<<"\n";
-
- // save point corresponding to the kink in the L curve as TGraph
- Double_t t[1],x[1],y[1];
- logTauX->GetKnot(iBest,t[0],x[0]);
- logTauY->GetKnot(iBest,t[0],y[0]);
- TGraph *bestLcurve=new TGraph(1,x,y);
- TGraph *bestLogTauX=new TGraph(1,t,x);
-
- //============================================================
- // extract unfolding results into histograms
-
- // set up a bin map, excluding underflow and overflow bins
- // the binMap relates the the output of the unfolding to the final
- // histogram bins
- Int_t *binMap=new Int_t[nGen+2];
- for(Int_t i=1;i<=nGen;i++) binMap[i]=i;
- binMap[0]=-1;
- binMap[nGen+1]=-1;
-
- TH1D *histMunfold=new TH1D("Unfolded",";mass(gen)",nGen,xminGen,xmaxGen);
- unfold.GetOutput(histMunfold,binMap);
- TH1D *histMdetFold=new TH1D("FoldedBack","mass(det)",nDet,xminDet,xmaxDet);
- unfold.GetFoldedOutput(histMdetFold);
-
- // store global correlation coefficients
- TH1D *histRhoi=new TH1D("rho_I","mass",nGen,xminGen,xmaxGen);
- unfold.GetRhoI(histRhoi,binMap);
-
- delete[] binMap;
- binMap=0;
-
- //=====================================================================
- // plot some histograms
- TCanvas *output = new TCanvas();
-
- // 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("*");
-
- return 0;
-}
diff --git a/tutorials/math/testUnfold3.C b/tutorials/math/testUnfold3.C
deleted file mode 100644
index 724ddb9dce7f721d0ee31e52519c95c18f5c72cb..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold3.C
+++ /dev/null
@@ -1,610 +0,0 @@
-/// \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
-///
-/// History:
-/// - Version 17.0, changeto 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
-///
-/// \macro_image
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#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: retreive unfolding results
-
- // get unfolding output
- // includes the statistical and background errors
- // but not the other systematic uncertainties
- TH1 *histUnfoldOutput=unfold.GetOutput("PT(unfold,stat+bgrerr)");
-
- // retreive error matrix of statistical errors
- TH2 *histEmatStat=unfold.GetEmatrixInput("unfolding stat error matrix");
-
- // retreive full error matrix
- // This includes all systematic errors
- TH2 *histEmatTotal=unfold.GetEmatrixTotal("unfolding total error matrix");
-
- // create two copies of the unfolded data, one with statistical errors
- // the other with total errors
- TH1 *histUnfoldStat=new TH1D("PT(unfold,staterr)",";Pt(gen)",
- nGen,xminGen,xmaxGen);
- TH1 *histUnfoldTotal=new TH1D("PT(unfold,totalerr)",";Pt(gen)",
- nGen,xminGen,xmaxGen);
- for(Int_t i=0;i<nGen+2;i++) {
- Double_t c=histUnfoldOutput->GetBinContent(i);
-
- // histogram with unfolded data and stat errors
- histUnfoldStat->SetBinContent(i,c);
- histUnfoldStat->SetBinError
- (i,TMath::Sqrt(histEmatStat->GetBinContent(i,i)));
-
- // histogram with unfolded data and total errors
- histUnfoldTotal->SetBinContent(i,c);
- histUnfoldTotal->SetBinError
- (i,TMath::Sqrt(histEmatTotal->GetBinContent(i,i)));
- }
-
- // create histogram with correlation matrix
- TH2D *histCorr=new TH2D("Corr(total)",";Pt(gen);Pt(gen)",
- nGen,xminGen,xmaxGen,nGen,xminGen,xmaxGen);
- for(Int_t i=0;i<nGen+2;i++) {
- Double_t ei,ej;
- ei=TMath::Sqrt(histEmatTotal->GetBinContent(i,i));
- if(ei<=0.0) continue;
- for(Int_t j=0;j<nGen+2;j++) {
- ej=TMath::Sqrt(histEmatTotal->GetBinContent(j,j));
- if(ej<=0.0) continue;
- histCorr->SetBinContent(i,j,histEmatTotal->GetBinContent(i,j)/ei/ej);
- }
- }
-
- // retreive bgr source 1
- TH1 *histDetNormBgr1=unfold.GetBackground("bgr1 normalized",
- "background1");
- TH1 *histDetNormBgrTotal=unfold.GetBackground("bgr total normalized");
-
- //========================
- // Step 7: plots
-
- TCanvas *output = new TCanvas();
- 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");
-
- //============================================================
- // 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/math/testUnfold4.C b/tutorials/math/testUnfold4.C
deleted file mode 100644
index 8bc350999557dbe928780665512b2a9b4656d514..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold4.C
+++ /dev/null
@@ -1,205 +0,0 @@
-/// \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.
-///
-/// History:
-/// - 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
-///
-/// \macro_image
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#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
-
- // switch off info messages
-#ifndef DEBUG
- gErrorIgnoreLevel = 1001;
-#endif
-
- 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 = new TCanvas();
- 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();
- }
-}
diff --git a/tutorials/math/testUnfold5a.C b/tutorials/math/testUnfold5a.C
deleted file mode 100644
index 7c986e601bcdc34948dc539fd1621c259818a79b..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold5a.C
+++ /dev/null
@@ -1,319 +0,0 @@
-/// \defgroup tutorial_unfold5 TUnfoldDensity and TUnfoldBinning test suite
-/// \ingroup tutorial_unfold
-///
-/// 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
-///
-/// \file
-/// \ingroup tutorial_unfold5
-/// \notebook -nodraw
-/// Version 17.0 example for multi-dimensional unfolding
-///
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#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(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/math/testUnfold5b.C b/tutorials/math/testUnfold5b.C
deleted file mode 100644
index 7ebe37d26af2aaec7f5ab2f3f7218a76ffbfa177..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold5b.C
+++ /dev/null
@@ -1,114 +0,0 @@
-/// \file
-/// \ingroup tutorial_unfold5
-/// \notebook -nodraw
-///
-/// Version 17.0 example for multi-dimensional unfolding
-///
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#include <iostream>
-#include <fstream>
-#include <TFile.h>
-#include "TUnfoldBinning.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
- );
- // 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();
-
- delete binningSchemes;
-}
diff --git a/tutorials/math/testUnfold5c.C b/tutorials/math/testUnfold5c.C
deleted file mode 100644
index f15f57d6aa9205f1e0d59faf25f41359571a6b98..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold5c.C
+++ /dev/null
@@ -1,211 +0,0 @@
-/// \file
-/// \ingroup tutorial_unfold5
-/// \notebook -nodraw
-/// Version 17.0 example for multi-dimensional unfolding
-///
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-
-#include <iostream>
-#include <map>
-#include <cmath>
-#include <TMath.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <TH1.h>
-#include "TUnfoldBinning.h"
-
-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 file
- // and save them to output file
-
- TFile *binningSchemes=new TFile("testUnfold5_binning.root");
-
- TUnfoldBinning *detectorBinning,*generatorBinning;
-
- outputFile->cd();
-
- binningSchemes->GetObject("detector",detectorBinning);
- binningSchemes->GetObject("generator",generatorBinning);
-
- delete binningSchemes;
-
- 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 bining scheme
- const TUnfoldBinning *detectordistribution=
- detectorBinning->FindNode("detectordistribution");
-
- const TUnfoldBinning *signalBinning=
- generatorBinning->FindNode("signal");
-
- const TUnfoldBinning *bgrBinning=
- generatorBinning->FindNode("background");
-
- // write binnig 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/math/testUnfold5d.C b/tutorials/math/testUnfold5d.C
deleted file mode 100644
index 1758e2aa619a0c30c785d7cd335346d81ae6c375..0000000000000000000000000000000000000000
--- a/tutorials/math/testUnfold5d.C
+++ /dev/null
@@ -1,216 +0,0 @@
-/// \file
-/// \ingroup tutorial_unfold5
-/// \notebook
-/// Version 17.0 example for multi-dimensional unfolding
-///
-/// \macro_image
-/// \macro_output
-/// \macro_code
-///
-/// \author Stefan Schmitt, DESY
-#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
-
-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);
-
- 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)
- unfold.SetInput(histDataReco /* ,0.0,1.0 */);
-
- // 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: retreive 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 = new TCanvas();
- 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]");
-
-}