-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathgenerateToyMC.C
119 lines (87 loc) · 3.46 KB
/
generateToyMC.C
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
/*
* generateToyMC.C
*
* Created on: 14 lis 2016
* Author: Leszek Kosarzewski
*/
#include <iostream>
#include "TString.h"
#include "TFile.h"
#include "TH1.h"
#include "TCanvas.h"
#include "TStyle.h"
#include "TROOT.h"
#include "style.C"
#ifndef __CINT__
#include "RooGlobalFunc.h"
#endif
#include "RooRealVar.h"
#include "RooDataHist.h"
#include "RooPlot.h"
#include "RooAddPdf.h"
#include "RooGaussian.h"
#include "RooExponential.h"
#include "RooHist.h"
using namespace std;
using namespace RooFit;
/////////////////////////////////////////////////////////////
void generateToyMC(TString fileName = "toyMC.root") {
style();
gStyle->SetOptFit(1);
gStyle->SetStatY(0.9);
gROOT->ForceStyle();
Double_t hist_minX = 5.0;
Double_t hist_maxX = 15.0;
TFile *file = new TFile(fileName, "RECREATE");
TH1D *histSB;
TH1D *histBkg;
TCanvas *cnv = new TCanvas("cnv", "cnv");
// Create parameters, set initial values and ranges
RooRealVar Bkg_A("Bkg_A","A", -0.5);
RooRealVar Bkg_B("Bkg_B","B", -0.05);
RooRealVar Sig_mu("Sig_mu","#mu", 8.0);
RooRealVar Sig_sigma("Sig_sigma","#sigma", 0.3);
RooRealVar Sig_mu_2("Sig_mu_2","#mu", 9.0);
RooRealVar Sig_sigma_2("Sig_sigma_2","#sigma", 0.4);
RooRealVar Sig_mu_3("Sig_mu_3","#mu", 12.0);
RooRealVar Sig_sigma_3("Sig_sigma_3","#sigma", 0.2);
RooRealVar Bkg_frac("Bkg_frac","Bkg_frac", 0.7);
RooRealVar Bkg_frac2("Bkg_frac","Bkg_frac", 0.2);
RooRealVar Sig_frac_1("Sig_frac_1","Sig_frac_1", 0.5);
RooRealVar Sig_frac_2("Sig_frac_2","Sig_frac_2", 0.6);
// Create variable, set range and it`s unit
RooRealVar mass("m_{ee}", "m_{ee}", hist_minX, hist_maxX, "GeV/c^{2}");
mass.setBins(100); // set number of bins for histogram
// Create PDFs
RooExponential *exp = new RooExponential("exp","Exp background", mass, Bkg_A);
RooExponential *exp2 = new RooExponential("exp2","Exp background 2", mass, Bkg_B);
RooGaussian *gauss = new RooGaussian("gauss", "Gaussian signal", mass, Sig_mu, Sig_sigma);
RooGaussian *gauss2 = new RooGaussian("gauss2", "Gaussian signal 2", mass, Sig_mu_2, Sig_sigma_2);
RooGaussian *gauss3 = new RooGaussian("gauss3", "Gaussian signal 3", mass, Sig_mu_3, Sig_sigma_3);
// Add PDFs together with recursive fractions (forces all ratios to be in [0,1])
//RooAddPdf *total = new RooAddPdf("total", "signal+background", RooArgList(*exp, *exp2, *gauss, *gauss2, *gauss3), RooArgList(Bkg_frac, Bkg_frac2, Sig_frac_1, Sig_frac_2), kTRUE);
RooAddPdf *total = new RooAddPdf("total", "signal+background", RooArgList(*exp, *gauss, *gauss2, *gauss3), RooArgList(Bkg_frac, Sig_frac_1, Sig_frac_2), kTRUE);
// Generate binned MC data - 100k events
RooDataHist *RooHistSig = total->generateBinned(mass, 100000);
RooDataHist *RooHistBkg = exp2->generateBinned(mass, 40000);
RooHistSig->add(*RooHistBkg);
// Export TH1D with 100 bins
histSB = (TH1D *)RooHistSig->createHistogram("histSB", mass, 100);
histBkg = (TH1D *)RooHistBkg->createHistogram("histBkg", mass, 100);
histSB->SetTitle("Signal; m_{ee} [GeV/c^2]; #frac{dN}{dm_{ee}} [#frac{c^2}{GeV}]");
histBkg->SetTitle("Background exp2; m_{ee} [GeV/c^2]; #frac{dN}{dm_{ee}} [#frac{c^2}{GeV}]");
// Plot simulated data
cnv->cd();
RooPlot* frame = mass.frame();
frame->SetTitle("Data");
RooHistSig->plotOn(frame);
RooHistBkg->plotOn(frame, MarkerColor(kRed));
frame->Draw();
//--------------------------------------------
file->cd();
histSB->Write("hSB");
histBkg->Write("hBkg");
cnv->cd()->SaveAs("toyMC.png");
cout<<"Finished!"<<endl;
//-----------------
}