1 #ifndef __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__
2 #define __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__
28 namespace JSUPERNOVA {
42 int n = in->GetNbinsX();
46 for (
int i = 0; i <
n; i++) {
47 out[i] = in->GetBinContent(i + 1);
89 TString rt_tag =
".RT_DET_SUM";
90 TString nc_tag =
".NC_DET_SUM";
91 TString rt2d_tag =
".RT2D_DET";
95 TIter iter(in->GetListOfKeys());
97 for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
99 const TString tag(key->GetName());
102 if (tag.EndsWith(rt_tag)) {
104 TString rt_hist_tag = tag;
105 TString run_tag = TString(tag, tag.Length() - rt_tag.Length());
107 TString nc_hist_tag = run_tag + nc_tag;
109 TH1D* RT = (TH1D*) in->Get(rt_hist_tag);
110 TH1D* NC = (TH1D*) in->Get(nc_hist_tag);
115 int n = nc_buf.size();
122 for (
int i = 0; i <
n; i++) {
125 nc.back().push_back( nc_buf[i] );
126 for (
int j = 100 * i;
j < 100 * (i + 1);
j++) {
127 rt.back().push_back( rt_buf[
j] );
133 TString rt2d_hist_tag = run_tag + rt2d_tag;
134 h2d_t* RT2D = (
h2d_t*) in->Get(rt2d_hist_tag);
137 RT2D->SetDirectory(0);
138 rt2d.push_back(RT2D);
179 void configureRatio(
const int inputNumberOfLines,
const int outputNumberOfLines) {
180 outputSize = inputNumberOfLines;
182 genRatio = outputNumberOfLines / inputNumberOfLines;
184 genScale = outputNumberOfLines / (1.0 * inputNumberOfLines * genRatio);
198 int k = rnd->Integer(rt.size());
200 int limit = rt[
k].size() - (genRatio * TWindow_ms);
202 int start = rnd->Integer(limit);
204 for (
int i = 0; i < genRatio; i++) {
206 for (
int j = 0;
j < TWindow_ms;
j++) {
208 int bin = start + (i * TWindow_ms) +
j;
210 out[0][
j] += genScale * rt[
k].at(bin );
211 out[1][
j] += genScale * nc[
k].at(bin / 100);
227 double mu = outputSize * 18 * domRate_Hz * 1.0e-3;
229 for (
unsigned i = 0; i < out[0].size(); i++) {
231 double count = rnd->Poisson(mu);
234 out[1][i] = outputSize * 18 * 31;
252 int nRuns = rt.size();
254 int k = rnd->Integer(nRuns);
259 int offset = rnd->Integer(100);
261 for (
int i = 0; i < genRatio; i++) {
263 int start = 100 * rnd->Integer(nc[k].size() - ceil(TWindow_ms / 100.0)) + offset;
265 for (
int j = 0;
j < TWindow_ms;
j++) {
269 out[0][
j] += genScale * rt[
k].at(bin );
270 out[1][
j] += genScale * nc[
k].at(bin / 100);
275 k = rnd->Integer(nRuns);
291 int ny = proto->GetNbinsY();
292 int ymin = proto->GetYaxis()->GetXmin();
293 int ymax = proto->GetYaxis()->GetXmax();
295 h2d_t* out =
new h2d_t(
"sample", NULL, TWindow_ms, 0, TWindow_ms, ny, ymin, ymax);
297 for (
int i = 0; i < genRatio; i++) {
299 for (
int j = 0;
j < TWindow_ms;
j++) {
301 int xbin = 1 + start + (i * TWindow_ms) +
j;
303 for (
int ybin = 1; ybin <= rt2d[
k]->GetNbinsY(); ybin++) {
305 double val = rt2d[
k]->GetBinContent(xbin, ybin);
307 double bcx = ((TAxis*)out->GetXaxis())->GetBinCenter(
j + 1);
308 double bcy = ((TAxis*)out->GetYaxis())->GetBinCenter(ybin);
310 out->Fill(bcx, bcy, val);
330 for (
int i = 0; i < genRatio; i++) {
332 for (
int j = 0;
j < TWindow_ms;
j++) {
334 int bin = start + (i * TWindow_ms) +
j;
336 out[
j] += genScale * nc[
k].at(bin / 100);
350 const int n = in.size();
354 for (
int i = 0; i <
n; i++) {
369 if (rt2d.size() > 0) {
371 int k = rnd->Integer(rt2d.size());
373 int limit = rt2d[
k]->GetNbinsX() - (genRatio * TWindow_ms);
375 int start = rnd->Integer(limit);
377 out = build_H2D(k, start);
381 if (rb != 1) { out->RebinX(rb); }
393 int n = in->GetNbinsX();
401 for (
int bin = 1; bin <=
n; bin++) {
403 double sg = 0, bg = 0, er = 0;
405 TH1D* h = in->ProjectionY(
"timeBin", bin, bin);
407 if (h->GetEntries() > 0) {
409 TF1
f(
"f",
"[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
411 f.SetParameter(0, h->GetMaximum());
412 f.SetParameter(1, h->GetMean());
413 f.SetParameter(2, h->GetRMS() * 0.25);
414 f.SetParameter(3, h->GetMinimum());
416 h->Fit(&f,
"Q",
"same");
418 bg = f.GetParameter(3) * h->GetNbinsX();
419 sg = h->GetSumOfWeights() - bg;
420 er = f.GetParError(3) * h->GetNbinsX();
447 int k = rnd->Integer(rt.size());
449 int limit = rt[
k].size() - (genRatio * TWindow_ms);
451 int start = rnd->Integer(limit);
453 h2d_t* sample = build_H2D(k, start);
455 if (rb != 1) { sample->RebinX(rb); }
459 out.push_back(fit[0]);
464 out.push_back(build_NC(k, start));
466 out.push_back(rebin(build_NC(k, start), rb));
469 out.push_back(fit[1]);
470 out.push_back(fit[2]);
*fatal Wrong number of arguments esac JCookie sh JRuns D $DETECTOR d sort n k
JL0BackgroundSimulator(TFile *in, bool twoDim=false)
Default constructor.
bg_type generate_fitted(int rb=1)
Generate fitted L1 sample.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
vector< vector< double > > bg_type
h2d_t * build_H2D(const int k, const int start)
Generate 2D histogram.
~JL0BackgroundSimulator()
Destructor.
bg_type generate_shuffled(bool randomizeRun=false)
Generate sample of L0 background The sampling of the L0 data is not sequential but random within the ...
vector< vector< double > > rt
static bg_type fit_H2D(h2d_t *in)
Fit 2D sample.
Class to emulate L0 background for an arbitrarily sized detector.
vector< vector< double > > nc
void setSeed(const UInt_t uSeed=0)
Set TRandom seed.
vector< double > build_NC(const int k, const int start)
Build NC sequence.
void configureRatio(const int inputNumberOfLines, const int outputNumberOfLines)
Configure generation ratio.
static vector< double > rebin(const vector< double > &in, const int rb)
Rebin vector.
alias put_queue eval echo n
void configureTimeWindow(const int T_ms)
Configure the duration of an output sample.
bg_type generate_poisson(const int domRate_Hz=500)
Generate pure poissionian background.
bg_type generate()
Generate sample of L0 background L0 data are randomly sampled from a single L0 dataset.
vector< double > loadHistogram(TH1D *in)
Load histogram values to vector, each bin is converted to an element.
h2d_t * generate_H2D(int rb=1)
Generate 2D sample.