arXiv:1611.01927 12th Conference on Quark Confinement and the Hadron Spectrum (Confinement XII)
This is an example of unfolding a two-dimensional distribution also using an auxiliary measurement to constrain some background
#include <iostream>
#include <map>
#include <cmath>
#define MASS1 0.511E-3
using namespace std;
class ToyEvent7 {
public:
void GenerateDataEvent(
TRandom *rnd);
void GenerateSignalEvent(
TRandom *rnd);
void GenerateBgrEvent(
TRandom *rnd);
inline Double_t GetMRec(
int i)
const {
return fMRec[i]; }
inline Double_t GetPtRec(
int i)
const {
return fPtRec[i]; }
inline Double_t GetEtaRec(
int i)
const {
return fEtaRec[i]; }
inline Double_t GetDiscriminator(
void)
const {
return fDiscriminator; }
inline Double_t GetPhiRec(
int i)
const {
return fPhiRec[i]; }
inline Bool_t IsTriggered(
void)
const {
return fIsTriggered; }
if(IsSignal()) return fMGen[i];
else return -1.0;
}
if(IsSignal()) return fPtGen[i];
else return -1.0;
}
inline Double_t GetEtaGen(
int i)
const {
if(IsSignal()) return fEtaGen[i];
else return 999.0;
}
inline Double_t GetPhiGen(
int i)
const {
if(IsSignal()) return fPhiGen[i];
else return 999.0;
}
inline Bool_t IsSignal(
void)
const {
return fIsSignal; }
protected:
public:
};
void testUnfold7a()
{
Int_t neventSigmc = 250000;
Int_t neventBgrmc = 100000;
Float_t etaRec[3],ptRec[3],phiRec[3],mRec[3],discr;
Float_t etaGen[3],ptGen[3],phiGen[3],mGen[3];
Int_t istriggered,issignal;
TFile *dataFile=
new TFile(
"testUnfold7_data.root",
"recreate");
dataTree->
Branch(
"etarec",etaRec,
"etarec[3]/F");
dataTree->
Branch(
"ptrec",ptRec,
"ptrec[3]/F");
dataTree->
Branch(
"phirec",phiRec,
"phirec[3]/F");
dataTree->
Branch(
"mrec",mRec,
"mrec[3]/F");
dataTree->
Branch(
"discr",&discr,
"discr/F");
dataTree->
Branch(
"istriggered",&istriggered,
"istriggered/I");
dataTree->
Branch(
"etagen",etaGen,
"etagen[3]/F");
dataTree->
Branch(
"ptgen",ptGen,
"ptgen[3]/F");
dataTree->
Branch(
"phigen",phiGen,
"phigen[3]/F");
dataTree->
Branch(
"mgen",mGen,
"mgen[3]/F");
dataTree->
Branch(
"issignal",&issignal,
"issignal/I");
cout<<"fill data tree\n";
for(int ievent=0;ievent<neventData;ievent++) {
event.GenerateDataEvent(g_rnd);
for(int i=0;i<3;i++) {
etaRec[i]=event.GetEtaRec(i);
ptRec[i]=event.GetPtRec(i);
phiRec[i]=event.GetPhiRec(i);
mRec[i]=event.GetMRec(i);
etaGen[i]=event.GetEtaGen(i);
ptGen[i]=event.GetPtGen(i);
phiGen[i]=event.GetPhiGen(i);
mGen[i]=event.GetMGen(i);
}
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
issignal=event.IsSignal() ? 1 : 0;
if(!(ievent%100000)) cout<<" data event "<<ievent<<"\n";
}
delete dataTree;
delete dataFile;
TFile *signalFile=
new TFile(
"testUnfold7_signal.root",
"recreate");
signalTree->
Branch(
"etarec",etaRec,
"etarec[3]/F");
signalTree->
Branch(
"ptrec",ptRec,
"ptrec[3]/F");
signalTree->
Branch(
"phirec",ptRec,
"phirec[3]/F");
signalTree->
Branch(
"mrec",mRec,
"mrec[3]/F");
signalTree->
Branch(
"discr",&discr,
"discr/F");
signalTree->
Branch(
"istriggered",&istriggered,
"istriggered/I");
signalTree->
Branch(
"etagen",etaGen,
"etagen[3]/F");
signalTree->
Branch(
"ptgen",ptGen,
"ptgen[3]/F");
signalTree->
Branch(
"phigen",phiGen,
"phigen[3]/F");
signalTree->
Branch(
"weight",&weight,
"weight/F");
signalTree->
Branch(
"mgen",mGen,
"mgen[3]/F");
cout<<"fill signal tree\n";
weight=ToyEvent7::kMCSignalFraction*muData0/neventSigmc;
for(int ievent=0;ievent<neventSigmc;ievent++) {
event.GenerateSignalEvent(g_rnd);
for(int i=0;i<3;i++) {
etaRec[i]=event.GetEtaRec(i);
ptRec[i]=event.GetPtRec(i);
phiRec[i]=event.GetPhiRec(i);
mRec[i]=event.GetMRec(i);
etaGen[i]=event.GetEtaGen(i);
ptGen[i]=event.GetPtGen(i);
phiGen[i]=event.GetPhiGen(i);
mGen[i]=event.GetMGen(i);
}
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
if(!(ievent%100000)) cout<<" signal event "<<ievent<<"\n";
}
delete signalTree;
delete signalFile;
TFile *bgrFile=
new TFile(
"testUnfold7_background.root",
"recreate");
bgrTree->
Branch(
"etarec",&etaRec,
"etarec[3]/F");
bgrTree->
Branch(
"ptrec",&ptRec,
"ptrec[3]/F");
bgrTree->
Branch(
"phirec",&phiRec,
"phirec[3]/F");
bgrTree->
Branch(
"mrec",&mRec,
"mrec[3]/F");
bgrTree->
Branch(
"discr",&discr,
"discr/F");
bgrTree->
Branch(
"istriggered",&istriggered,
"istriggered/I");
bgrTree->
Branch(
"weight",&weight,
"weight/F");
cout<<"fill background tree\n";
weight=(1.-ToyEvent7::kMCSignalFraction)*muData0/neventBgrmc;
for(int ievent=0;ievent<neventBgrmc;ievent++) {
event.GenerateBgrEvent(g_rnd);
for(int i=0;i<3;i++) {
etaRec[i]=event.GetEtaRec(i);
ptRec[i]=event.GetPtRec(i);
phiRec[i]=event.GetPhiRec(i);
}
discr=event.GetDiscriminator();
istriggered=event.IsTriggered() ? 1 : 0;
if(!(ievent%100000)) cout<<" background event "<<ievent<<"\n";
}
delete bgrTree;
delete bgrFile;
}
Double_t ToyEvent7::kDataSignalFraction=0.75;
Double_t ToyEvent7::kMCSignalFraction=0.75;
void ToyEvent7::GenerateDataEvent(
TRandom *rnd) {
fIsSignal=rnd->
Uniform()<kDataSignalFraction;
if(IsSignal()) {
GenerateSignalKinematics(rnd,
kTRUE);
} else {
GenerateBgrKinematics(rnd,
kTRUE);
}
GenerateReco(rnd);
}
void ToyEvent7::GenerateSignalEvent(
TRandom *rnd) {
fIsSignal=true;
GenerateSignalKinematics(rnd,
kFALSE);
GenerateReco(rnd);
}
void ToyEvent7::GenerateBgrEvent(
TRandom *rnd) {
fIsSignal=false;
GenerateBgrKinematics(rnd,
kFALSE);
GenerateReco(rnd);
}
void ToyEvent7::GenerateSignalKinematics(
TRandom *rnd,
Bool_t isData) {
double M0=91.1876;
do {
} while(fMGen[2]<=0.0);
double N_ETA=3.0;
double MU_PT=5.;
double SIGMA_PT=2.0;
double DECAY_A=0.2;
if(isData) {
MU_PT=6.;
SIGMA_PT=1.8;
}
if(rnd->
Uniform(-1.,1.)<0.) fEtaGen[2] *= -1.;
do {
fPtGen[2]=rnd->
Landau(MU_PT,SIGMA_PT);
} while((fPtGen[2]<=0.0)||(fPtGen[2]>500.));
double costh;
do {
costh=
r*(1.+DECAY_A*
r*
r);
}
while(
fabs(costh)>=1.0);
double pz=ptot*costh;
p[0].SetXYZT(px,py,pz,
e);
p[1].SetXYZT(-px,-py,-pz,
e);
for(int i=0;i<2;i++) {
}
for(int i=0;i<3;i++) {
}
}
void ToyEvent7::GenerateBgrKinematics(
TRandom *rnd,
Bool_t isData) {
for(int i=0;i<3;i++) {
fPtGen[i]=0.0;
fEtaGen[i]=0.0;
fPhiGen[i]=0.0;
}
for(int i=0;i<2;i++) {
}
for(int i=0;i<3;i++) {
}
}
void ToyEvent7::GenerateReco(
TRandom *rnd) {
if(fIsSignal) {
for(int i=0;i<2;i++) {
+1.0*coshEta
+0.01*eGen;
do {
eRec=rnd->
Gaus(eGen,sigmaE);
} while(eRec<=0.0);
p[i].SetPtEtaPhiM(eRec/(expEta+1./expEta),
rnd->
Gaus(fEtaGen[i],sigmaEta),
MASS1);
}
for(int i=0;i<3;i++) {
}
}
if(fIsSignal) {
do {
Double_t tauDiscr=0.08-0.04/(1.+fPtRec[2]/10.0);
fDiscriminator=1.0-rnd->
Exp(tauDiscr)+rnd->
Gaus(0.,sigmaDiscr);
} while((fDiscriminator<=0.)||(fDiscriminator>=1.));
} else {
do {
Double_t tauDiscr=0.15-0.05/(1.+fPtRec[2]/5.0)+0.1*fEtaRec[2];
Double_t sigmaDiscr=0.02+0.01*fEtaRec[2];
fDiscriminator=rnd->
Exp(tauDiscr)+rnd->
Gaus(0.,sigmaDiscr);
} while((fDiscriminator<=0.)||(fDiscriminator>=1.));
}
fIsTriggered=false;
for(int i=0;i<2;i++) {
}
}
winID h TVirtualViewer3D TVirtualGLPainter p
Option_t Option_t TPoint TPoint const char GetTextMagnitude GetFillStyle GetLineColor GetLineWidth GetMarkerStyle GetTextAlign GetTextColor GetTextSize void char Point_t Rectangle_t WindowAttributes_t Float_t r
A ROOT file is composed of a header, followed by consecutive data records (TKey instances) with a wel...
void SetPtEtaPhiM(Double_t pt, Double_t eta, Double_t phi, Double_t m)
Random number generator class based on M.
This is the base class for the ROOT Random number generators.
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...
virtual Int_t Poisson(Double_t mean)
Generates a random integer N according to a Poisson law.
virtual Double_t Exp(Double_t tau)
Returns an exponential deviate.
virtual Double_t Uniform(Double_t x1=1)
Returns a uniform deviate on the interval (0, x1).
virtual Double_t Landau(Double_t mean=0, Double_t sigma=1)
Generate a random number following a Landau distribution with location parameter mu and scale paramet...
virtual Double_t BreitWigner(Double_t mean=0, Double_t gamma=1)
Return a number distributed following a BreitWigner function with mean and gamma.
A TTree represents a columnar dataset.
virtual Int_t Fill()
Fill all branches.
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.
Int_t Write(const char *name=nullptr, Int_t option=0, Int_t bufsize=0) override
Write this object to the current directory.
RVec< PromoteTypes< T0, T1 > > remainder(const T0 &x, const RVec< T1 > &v)
RVec< PromoteType< T > > cos(const RVec< T > &v)
RVec< PromoteType< T > > sin(const RVec< T > &v)
LVector boost(const LVector &v, const BoostVector &b)
Boost a generic Lorentz Vector class using a generic 3D Vector class describing the boost The only re...
VecExpr< UnaryOp< Fabs< T >, VecExpr< A, T, D >, T >, T, D > fabs(const VecExpr< A, T, D > &rhs)
Double_t Exp(Double_t x)
Returns the base-e exponential function of x, which is e raised to the power x.
Double_t Sqrt(Double_t x)
Returns the square root of x.
LongDouble_t Power(LongDouble_t x, LongDouble_t y)
Returns x raised to the power y.
Short_t Abs(Short_t d)
Returns the absolute value of parameter Short_t d.
static uint64_t sum(uint64_t i)