Logo ROOT  

Reference Guide
testUnfold5a.C File Reference
Tutorials » TUnfold tutorials

Detailed Description

View in nbviewer Open in SWAN Test program for the classes TUnfoldDensity and TUnfoldBinning

A toy test of the TUnfold package

This is an example of unfolding a two-dimensional distribution also using an auxiliary measurement to constrain some background

The example comprises several macros

  • testUnfold5a.C create root files with TTree objects for signal, background and data
    • write files testUnfold5_signal.root testUnfold5_background.root testUnfold5_data.root
  • testUnfold5b.C create a root file with the TUnfoldBinning objects
    • write file testUnfold5_binning.root
  • testUnfold5c.C loop over trees and fill histograms based on the TUnfoldBinning objects
    • read testUnfold5_binning.root testUnfold5_signal.root testUnfold5_background.root testUnfold5_data.root
    • write testUnfold5_histograms.root
  • testUnfold5d.C run the unfolding
    • read testUnfold5_histograms.root
    • write testUnfold5_result.root testUnfold5_result.ps
fill data tree
data event 0
data event 100000
fill signal tree
signal event 0
signal event 100000
signal event 200000
signal event 300000
signal event 400000
signal event 500000
signal event 600000
signal event 700000
signal event 800000
signal event 900000
signal event 1000000
signal event 1100000
signal event 1200000
signal event 1300000
signal event 1400000
signal event 1500000
signal event 1600000
signal event 1700000
signal event 1800000
signal event 1900000
fill background tree
background event 0
background event 100000
background event 200000
background event 300000
background event 400000
background event 500000
background event 600000
background event 700000
background event 800000
background event 900000
background event 1000000
background event 1100000
background event 1200000
background event 1300000
background event 1400000
background event 1500000
background event 1600000
background event 1700000
background event 1800000
background event 1900000
#include <iostream>
#include <map>
#include <cmath>
#include <TMath.h>
#include <TRandom3.h>
#include <TFile.h>
#include <TTree.h>
using namespace std;
TRandom *g_rnd=0;
class ToyEvent {
public:
void GenerateDataEvent(TRandom *rnd);
void GenerateSignalEvent(TRandom *rnd);
void GenerateBgrEvent(TRandom *rnd);
// reconstructed quantities
inline Double_t GetPtRec(void) const { return fPtRec; }
inline Double_t GetEtaRec(void) const { return fEtaRec; }
inline Double_t GetDiscriminator(void) const {return fDiscriminator; }
inline Bool_t IsTriggered(void) const { return fIsTriggered; }
// generator level quantities
inline Double_t GetPtGen(void) const {
if(IsSignal()) return fPtGen;
else return -1.0;
}
inline Double_t GetEtaGen(void) const {
if(IsSignal()) return fEtaGen;
else return 999.0;
}
inline Bool_t IsSignal(void) const { return fIsSignal; }
protected:
void GenerateSignalKinematics(TRandom *rnd,Bool_t isData);
void GenerateBgrKinematics(TRandom *rnd,Bool_t isData);
void GenerateReco(TRandom *rnd);
// reconstructed quantities
Double_t fPtRec;
Double_t fEtaRec;
Double_t fDiscriminator;
Bool_t fIsTriggered;
// generated quantities
Double_t fPtGen;
Double_t fEtaGen;
Bool_t fIsSignal;
static Double_t kDataSignalFraction;
};
void testUnfold5a()
{
// random generator
g_rnd=new TRandom3();
// data and MC number of events
const Int_t neventData = 20000;
Double_t const neventSignalMC =2000000;
Double_t const neventBgrMC =2000000;
Float_t etaRec,ptRec,discr,etaGen,ptGen;
Int_t istriggered,issignal;
//==================================================================
// Step 1: generate data TTree
TFile *dataFile=new TFile("testUnfold5_data.root","recreate");
TTree *dataTree=new TTree("data","event");
dataTree->Branch("etarec",&etaRec,"etarec/F");
dataTree->Branch("ptrec",&ptRec,"ptrec/F");
dataTree->Branch("discr",&discr,"discr/F");
// for real data, only the triggered events are available
dataTree->Branch("istriggered",&istriggered,"istriggered/I");
// data truth parameters
dataTree->Branch("etagen",&etaGen,"etagen/F");
dataTree->Branch("ptgen",&ptGen,"ptgen/F");
dataTree->Branch("issignal",&issignal,"issignal/I");
cout<<"fill data tree\n";
Int_t nEvent=0,nTriggered=0;
while(nTriggered<neventData) {
ToyEvent event;
event.GenerateDataEvent(g_rnd);
etaRec=event.GetEtaRec();
ptRec=event.GetPtRec();
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
etaGen=event.GetEtaGen();
ptGen=event.GetPtGen();
issignal=event.IsSignal() ? 1 : 0;
dataTree->Fill();
if(!(nEvent%100000)) cout<<" data event "<<nEvent<<"\n";
if(istriggered) nTriggered++;
nEvent++;
}
dataTree->Write();
delete dataTree;
delete dataFile;
//==================================================================
// Step 2: generate signal TTree
TFile *signalFile=new TFile("testUnfold5_signal.root","recreate");
TTree *signalTree=new TTree("signal","event");
signalTree->Branch("etarec",&etaRec,"etarec/F");
signalTree->Branch("ptrec",&ptRec,"ptrec/F");
signalTree->Branch("discr",&discr,"discr/F");
signalTree->Branch("istriggered",&istriggered,"istriggered/I");
signalTree->Branch("etagen",&etaGen,"etagen/F");
signalTree->Branch("ptgen",&ptGen,"ptgen/F");
cout<<"fill signal tree\n";
for(int ievent=0;ievent<neventSignalMC;ievent++) {
ToyEvent event;
event.GenerateSignalEvent(g_rnd);
etaRec=event.GetEtaRec();
ptRec=event.GetPtRec();
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
etaGen=event.GetEtaGen();
ptGen=event.GetPtGen();
if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
signalTree->Fill();
}
signalTree->Write();
delete signalTree;
delete signalFile;
// ==============================================================
// Step 3: generate background MC TTree
TFile *bgrFile=new TFile("testUnfold5_background.root","recreate");
TTree *bgrTree=new TTree("background","event");
bgrTree->Branch("etarec",&etaRec,"etarec/F");
bgrTree->Branch("ptrec",&ptRec,"ptrec/F");
bgrTree->Branch("discr",&discr,"discr/F");
bgrTree->Branch("istriggered",&istriggered,"istriggered/I");
cout<<"fill background tree\n";
for(int ievent=0;ievent<neventBgrMC;ievent++) {
ToyEvent event;
event.GenerateBgrEvent(g_rnd);
etaRec=event.GetEtaRec();
ptRec=event.GetPtRec();
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
bgrTree->Fill();
}
bgrTree->Write();
delete bgrTree;
delete bgrFile;
}
Double_t ToyEvent::kDataSignalFraction=0.8;
void ToyEvent::GenerateDataEvent(TRandom *rnd) {
fIsSignal=rnd->Uniform()<kDataSignalFraction;
if(IsSignal()) {
GenerateSignalKinematics(rnd,kTRUE);
} else {
GenerateBgrKinematics(rnd,kTRUE);
}
GenerateReco(rnd);
}
void ToyEvent::GenerateSignalEvent(TRandom *rnd) {
fIsSignal=1;
GenerateSignalKinematics(rnd,kFALSE);
GenerateReco(rnd);
}
void ToyEvent::GenerateBgrEvent(TRandom *rnd) {
fIsSignal=0;
GenerateBgrKinematics(rnd,kFALSE);
GenerateReco(rnd);
}
void ToyEvent::GenerateSignalKinematics(TRandom *rnd,Bool_t isData) {
Double_t e_T0=0.5;
Double_t e_T0_eta=0.0;
Double_t e_n=2.0;
Double_t e_n_eta=0.0;
Double_t eta_p2=0.0;
Double_t etaMax=4.0;
if(isData) {
e_T0=0.6;
e_n=2.5;
e_T0_eta=0.05;
e_n_eta=-0.05;
eta_p2=1.5;
}
if(eta_p2>0.0) {
fEtaGen=TMath::Power(rnd->Uniform(),eta_p2)*etaMax;
if(rnd->Uniform()>=0.5) fEtaGen= -fEtaGen;
} else {
fEtaGen=rnd->Uniform(-etaMax,etaMax);
}
Double_t n=e_n + e_n_eta*fEtaGen;
Double_t T0=e_T0 + e_T0_eta*fEtaGen;
fPtGen=(TMath::Power(rnd->Uniform(),1./(1.-n))-1.)*T0;
/* static int print=100;
if(print) {
cout<<fEtaGen
<<" "<<fPtGen
<<"\n";
print--;
} */
}
void ToyEvent::GenerateBgrKinematics(TRandom *rnd,Bool_t isData) {
fPtGen=0.0;
fEtaGen=0.0;
fPtRec=rnd->Exp(isData ? 2.5 : 2.5);
fEtaRec=rnd->Uniform(-3.,3.);
}
void ToyEvent::GenerateReco(TRandom *rnd) {
if(fIsSignal) {
Double_t expEta=TMath::Exp(fEtaGen);
Double_t eGen=fPtGen*(expEta+1./expEta);
Double_t sigmaE=0.1*TMath::Sqrt(eGen)+(0.01+0.002*TMath::Abs(fEtaGen))
*eGen;
Double_t eRec;
do {
eRec=rnd->Gaus(eGen,sigmaE);
} while(eRec<=0.0);
Double_t sigmaEta=0.1+0.02*fEtaGen;
fEtaRec=rnd->Gaus(fEtaGen,sigmaEta);
fPtRec=eRec/(expEta+1./expEta);
do {
Double_t tauDiscr=0.08-0.04/(1.+fPtRec/10.0);
Double_t sigmaDiscr=0.01;
fDiscriminator=1.0-rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
} while((fDiscriminator<=0.)||(fDiscriminator>=1.));
/* static int print=100;
if(print) {
cout<<fEtaGen<<" "<<fPtGen
<<" -> "<<fEtaRec<<" "<<fPtRec
<<"\n";
print--;
} */
} else {
do {
Double_t tauDiscr=0.15-0.05/(1.+fPtRec/5.0)+0.1*fEtaRec;
Double_t sigmaDiscr=0.02+0.01*fEtaRec;
fDiscriminator=rnd->Exp(tauDiscr)+rnd->Gaus(0.,sigmaDiscr);
} while((fDiscriminator<=0.)||(fDiscriminator>=1.));
}
fIsTriggered=(rnd->Uniform()<1./(TMath::Exp(-fPtRec+3.5)+1.));
}
Int_t
int Int_t
Definition: RtypesCore.h:41
kFALSE
const Bool_t kFALSE
Definition: RtypesCore.h:88
Bool_t
bool Bool_t
Definition: RtypesCore.h:59
Double_t
double Double_t
Definition: RtypesCore.h:55
Float_t
float Float_t
Definition: RtypesCore.h:53
kTRUE
const Bool_t kTRUE
Definition: RtypesCore.h:87
TFile
A ROOT file is a suite of consecutive data records (TKey instances) with a well defined format.
Definition: TFile.h:48
TRandom3
Random number generator class based on M.
Definition: TRandom3.h:27
TRandom
This is the base class for the ROOT Random number generators.
Definition: TRandom.h:27
TRandom::Gaus
virtual Double_t Gaus(Double_t mean=0, Double_t sigma=1)
Samples a random number from the standard Normal (Gaussian) Distribution with the given mean and sigm...
Definition: TRandom.cxx:263
TRandom::Exp
virtual Double_t Exp(Double_t tau)
Returns an exponential deviate.
Definition: TRandom.cxx:240
TRandom::Uniform
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
Definition: TRandom.cxx:635
TTree
A TTree represents a columnar dataset.
Definition: TTree.h:72
TTree::Fill
virtual Int_t Fill()
Fill all branches.
Definition: TTree.cxx:4487
TTree::Branch
TBranch * Branch(const char *name, T *obj, Int_t bufsize=32000, Int_t splitlevel=99)
Add a new branch, and infer the data type from the type of obj being passed.
Definition: TTree.h:341
TTree::Write
virtual Int_t Write(const char *name=0, Int_t option=0, Int_t bufsize=0)
Write this object to the current directory.
Definition: TTree.cxx:9485
n
const Int_t n
Definition: legend1.C:16
ROOT::Math::etaMax
T etaMax()
Function providing the maximum possible value of pseudorapidity for a non-zero rho,...
Definition: etaMax.h:50
TMath::Exp
Double_t Exp(Double_t x)
Definition: TMath.h:717
TMath::Sqrt
Double_t Sqrt(Double_t x)
Definition: TMath.h:681
TMath::Power
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Definition: TMath.h:725
TMath::Abs
Short_t Abs(Short_t d)
Definition: TMathBase.h:120

Version 17.6, in parallel to changes in TUnfold

History:

  • Version 17.5, in parallel to changes in TUnfold
  • Version 17.4, in parallel to changes in TUnfold
  • Version 17.3, in parallel to changes in TUnfold
  • Version 17.2, in parallel to changes in TUnfold
  • Version 17.1, in parallel to changes in TUnfold
  • Version 17.0 example for multi-dimensional unfolding

This file is part of TUnfold.

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

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

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

Author
Stefan Schmitt DESY, 14.10.2008

Definition in file testUnfold5a.C.