Generate Monte Carlo and Data events The events consist of
background
The signal is a resonance. It is generated with a Breit-Wigner, smeared by a Gaussian
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
fit the unfolded distribution, including the correlation matrix
save six plots to a file testUnfold1.ps
Version 17.0, updated for using the classes TUnfoldDensity, TUnfoldBinning
History:
Author: Stefan Schmitt, DESY
This notebook tutorial was automatically generated with ROOTBOOK-izer (Beta) from the macro found in the ROOT repository on Thursday, January 19, 2017 at 04:34 PM.
In [ ]:
%%cpp -d
#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"
In [ ]:
using namespace std;
TRandom *rnd=0;
TH2 *gHistInvEMatrix;
TVirtualFitter *gFitter=0;
A helper function is created:
In [ ]:
%%cpp -d
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);
}
A helper function is created:
In [ ]:
%%cpp -d
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
In [ ]:
%%cpp -d
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
In [ ]:
%%cpp -d
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);
}
}
Switch on histogram errors
In [ ]:
TH1::SetDefaultSumw2();
Show fit result
In [ ]:
gStyle->SetOptFit(1111);
Random generator
In [ ]:
rnd=new TRandom3();
Data and mc luminosity, cross-section
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
Int_t nScan=30;
Use automatic l-curve scan: start with taumin=taumax=0.0
In [ ]:
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)
In [ ]:
#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
In [ ]:
iBest=unfold.ScanLcurve(nScan,tauMin,tauMax,&lCurve,&logTauX,&logTauY);
If required, switch to previous log-level
In [ ]:
#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
In [ ]:
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
In [ ]:
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
In [ ]:
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
In [ ]:
TH1 *histMunfold=unfold.GetOutput("Unfolded");
Get unfolding result, folded back
In [ ]:
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
In [ ]:
TH2 *histEmatTotal=unfold.GetEmatrixTotal("EmatTotal");
Create data histogram with the total errors
In [ ]:
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
In [ ]:
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
In [ ]:
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.
In [ ]:
histMunfold->Fit(bw,"UE");
plot some histograms
In [ ]:
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
In [ ]:
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)
In [ ]:
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)
In [ ]:
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
In [ ]:
output.cd(4);
histRhoi->Draw();
Show tau as a function of chi**2
In [ ]:
output.cd(5);
logTauX->Draw();
bestLogTauLogChi2->SetMarkerColor(kRed);
bestLogTauLogChi2->Draw("*");
Show the l curve
In [ ]:
output.cd(6);
lCurve->Draw("AL");
bestLcurve->SetMarkerColor(kRed);
bestLcurve->Draw("*");
output.SaveAs("testUnfold1.ps");
return 0;
Draw all canvases
In [ ]:
gROOT->GetListOfCanvases()->Draw()