Jpp
JL0BackgroundSimulator.hh
Go to the documentation of this file.
1 #ifndef __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__
2 #define __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__
3 
4 #include <vector>
5 
6 #include "TFile.h"
7 #include "TKey.h"
8 
9 #include "TH1D.h"
10 #include "TH2F.h"
11 #include "TF1.h"
12 
13 #include "TRandom3.h"
14 
15 #include "JRipple.hh"
16 
17 /**
18  * \file
19  *
20  * This file provides a class to emulate an L0 background on the ms time scale for an arbitrarily sized detector.
21  * The emulation is based on the sampling of real L0 data (in principle from a small-sized detector) through the JRipple application.
22  * Samples of given duration are taken in a contiguous sequence and summed bin-by-bin as they were coming from different DUs,
23  * in this way time-correlations are used to emulate space-correlations which are not possible to probe due to the detector limited size.
24  *
25  * \author mlincett
26  */
27 
28 namespace JSUPERNOVA {
29 
30  using namespace std;
31 
33 
34  /**
35  * Load histogram values to vector, each bin is converted to an element.
36  * No other information is preserved (binning, axis semantics, under/overflows).
37  *
38  * \param in input histogram
39  * \return vector with one value per bin
40  */
42  int n = in->GetNbinsX();
43 
44  vector<double> out(n);
45 
46  for (int i = 0; i < n; i++) {
47  out[i] = in->GetBinContent(i + 1);
48  }
49 
50  return out;
51  }
52 
53  /**
54  * Class to emulate L0 background for an arbitrarily sized detector.
55  */
56 
58  private:
61 
63 
65  int genRatio;
67  // int rebin;
68 
69  double genScale;
70 
71  TRandom* rnd;
72 
73  public:
74 
75  /**
76  * Default constructor
77  *
78  * The input file contains a sequence of histograms with sampled L0 data;
79  * Each histogram is loaded as a vector, which will be considered as an individual L0 dataset.
80  *
81  * \param in input file from JRipple
82  */
83  JL0BackgroundSimulator(TFile* in, bool twoDim = false) {
84 
85  TWindow_ms = 100;
86  genRatio = 115;
87  genScale = 1;
88  // rebin = 1;
89 
90  TString rt_tag = ".RT_DET_SUM";
91  TString nc_tag = ".NC_DET_SUM";
92  TString rt2d_tag = ".RT2D_DET";
93 
94  rnd = new TRandom();
95 
96  TIter iter(in->GetListOfKeys());
97 
98  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
99 
100  const TString tag(key->GetName());
101 
102  // match
103  if (tag.EndsWith(rt_tag)) {
104 
105  TString rt_hist_tag = tag;
106  TString run_tag = TString(tag, tag.Length() - rt_tag.Length());
107 
108  TString nc_hist_tag = run_tag + nc_tag;
109 
110  TH1D* RT = (TH1D*) in->Get(rt_hist_tag);
111  TH1D* NC = (TH1D*) in->Get(nc_hist_tag);
112 
113  vector<double> rt_buf = loadHistogram(RT);
114  vector<double> nc_buf = loadHistogram(NC);
115 
116  int n = nc_buf.size();
117 
118  rt.push_back( vector<double>() );
119  nc.push_back( vector<double>() );
120 
121  // build parallel vectors
122  // 1 NC bin <=> 100 RT bins
123  for (int i = 0; i < n; i++) {
124  // empty slices are discarded
125  if (nc_buf[i] > 0) {
126  nc.back().push_back( nc_buf[i] );
127  for (int j = 100 * i; j < 100 * (i + 1); j++) {
128  rt.back().push_back( rt_buf[j] );
129  }
130  }
131  }
132 
133  if (twoDim) {
134  TString rt2d_hist_tag = run_tag + rt2d_tag;
135  h2d_t* RT2D = (h2d_t*) in->Get(rt2d_hist_tag);
136 
137  if (RT2D != NULL) {
138  RT2D->SetDirectory(0);
139  rt2d.push_back(RT2D);
140  }
141  }
142 
143  }
144 
145  }
146 
147  }
148 
149  /**
150  * Destructor
151  */
153  delete rnd;
154  }
155 
156  /**
157  * Configure the duration of an output sample
158  *
159  * \param T_ms sample duration in ms
160  */
161  void configureTimeWindow(const int T_ms) {
162  TWindow_ms = T_ms;
163  }
164 
165  /**
166  * Set TRandom seed
167  *
168  * \param uSeed seed
169  */
170  void setSeed(const UInt_t uSeed = 0) {
171  rnd->SetSeed(uSeed);
172  }
173 
174  /**
175  * Set rebinning ratio
176  * \param rebin
177 
178  void setRebin(const int rb) {
179  rebin = rb;
180  }
181  */
182 
183  /**
184  * Configure generation ratio
185  *
186  * \param inputSize number of lines of input detector
187  * \param outputSize number of lines of output detector
188  */
189  void configureRatio(const int iSize, const int oSize) {
190  outputSize = oSize;
191  // ratio is integer by design
192  genRatio = oSize / iSize;
193  // scaling factor to compensate the ratio remainder
194  genScale = oSize / (1.0 * iSize * genRatio);
195  }
196 
197 
198  /**
199  * Generate sample of L0 background
200  * L0 data are randomly sampled from a single L0 dataset.
201  *
202  * \return 2D vector with hit count and number of active channels as a function of time
203  */
205 
206  bg_type out(2, vector<double>(TWindow_ms));
207 
208  int k = rnd->Integer(rt.size());
209 
210  int limit = rt[k].size() - (genRatio * TWindow_ms);
211 
212  int start = rnd->Integer(limit);
213 
214  for (int i = 0; i < genRatio; i++) {
215 
216  for (int j = 0; j < TWindow_ms; j++) {
217 
218  int bin = start + (i * TWindow_ms) + j;
219 
220  out[0][j] += genScale * rt[k].at(bin );
221  out[1][j] += genScale * nc[k].at(bin / 100);
222 
223  }
224 
225  }
226 
227  return out;
228  }
229 
230  /**
231  * Generate pure poissionian background
232  */
233  bg_type generate_poisson(const int domRate_Hz = 500) {
234 
235  bg_type out(2, vector<double>(TWindow_ms));
236 
237  double mu = outputSize * 18 * domRate_Hz * 1.0e-3;
238 
239  for (unsigned i = 0; i < out[0].size(); i++) {
240 
241  double count = rnd->Poisson(mu);
242 
243  out[0][i] = count;
244  out[1][i] = outputSize * 18 * 31;
245 
246  }
247 
248  return out;
249 
250  }
251 
252 
253  /**
254  * Generate sample of L0 background
255  * The sampling of the L0 data is not sequential but random within the L0 dataset
256  * \return 2D vector with hit count and number of active channels as a function of time
257  */
258  bg_type generate_shuffled(bool randomizeRun = false) {
259 
260  bg_type out(2, vector<double>(TWindow_ms));
261 
262  int nRuns = rt.size();
263 
264  int k = rnd->Integer(nRuns);
265 
266  // a series n = 'genRatio' elements is sampled randomly within the dataset
267  // all the elements start with the same offset to the timeslice beginning
268 
269  int offset = rnd->Integer(100);
270 
271  for (int i = 0; i < genRatio; i++) {
272 
273  int start = 100 * rnd->Integer(nc[k].size() - ceil(TWindow_ms / 100.0)) + offset;
274 
275  for (int j = 0; j < TWindow_ms; j++) {
276 
277  int bin = start + j;
278 
279  out[0][j] += genScale * rt[k].at(bin );
280  out[1][j] += genScale * nc[k].at(bin / 100);
281 
282  }
283 
284  if (randomizeRun) {
285  k = rnd->Integer(nRuns);
286  }
287 
288  }
289 
290  return out;
291  }
292 
293 
294  /**
295  * Generate 2D histogram
296  */
297  h2d_t* build_H2D (const int k, const int start) {
298 
299  h2d_t* proto = rt2d[k];
300 
301  int ny = proto->GetNbinsY();
302  int ymin = proto->GetYaxis()->GetXmin();
303  int ymax = proto->GetYaxis()->GetXmax();
304 
305  h2d_t* out = new h2d_t("sample", NULL, TWindow_ms, 0, TWindow_ms, ny, ymin, ymax);
306 
307  for (int i = 0; i < genRatio; i++) {
308 
309  for (int j = 0; j < TWindow_ms; j++) {
310 
311  int xbin = 1 + start + (i * TWindow_ms) + j;
312 
313  for (int ybin = 1; ybin <= rt2d[k]->GetNbinsY(); ybin++) {
314 
315  double val = rt2d[k]->GetBinContent(xbin, ybin);
316 
317  double bcx = ((TAxis*)out->GetXaxis())->GetBinCenter(j + 1);
318  double bcy = ((TAxis*)out->GetYaxis())->GetBinCenter(ybin);
319 
320  out->Fill(bcx, bcy, val);
321 
322  }
323 
324  }
325 
326  }
327 
328  return out;
329 
330  }
331 
332 
333  /**
334  * Build NC sequence
335  */
336  vector<double> build_NC (const int k, const int start) {
337 
338  vector<double> out(TWindow_ms);
339 
340  for (int i = 0; i < genRatio; i++) {
341 
342  for (int j = 0; j < TWindow_ms; j++) {
343 
344  int bin = start + (i * TWindow_ms) + j;
345 
346  out[j] += genScale * nc[k].at(bin / 100);
347 
348  }
349 
350  }
351 
352  return out;
353  }
354 
355  /**
356  * Rebin vector
357  */
358  static vector<double> rebin(const vector<double> &in, const int rb) {
359 
360  const int n = in.size();
361 
362  vector<double> out(n/rb);
363 
364  for (int i = 0; i < n; i++) {
365  out[i/rb] += in[i];
366  }
367 
368  return out;
369 
370  }
371 
372  /**
373  * Generate 2D sample
374  */
375  h2d_t* generate_H2D(int rb = 1) {
376 
377  h2d_t* out = NULL;
378 
379  if (rt2d.size() > 0) {
380 
381  int k = rnd->Integer(rt2d.size());
382 
383  int limit = rt2d[k]->GetNbinsX() - (genRatio * TWindow_ms);
384 
385  int start = rnd->Integer(limit);
386 
387  out = build_H2D(k, start);
388 
389  }
390 
391  if (rb != 1) { out->RebinX(rb); }
392 
393  return out;
394 
395  }
396 
397 
398  /**
399  * Fit 2D sample
400  */
401  static bg_type fit_H2D(h2d_t* in) {
402 
403  int n = in->GetNbinsX();
404 
405  // 0 -> signal
406  // 1 -> fitted background
407  // 2 -> background fit error
408 
409  bg_type out(3, vector<double>(n));
410 
411  for (int bin = 1; bin <= n; bin++) {
412 
413  double sg = 0, bg = 0, er = 0;
414 
415  TH1D* h = in->ProjectionY("timeBin", bin, bin);
416 
417  if (h->GetEntries() > 0) {
418 
419  TF1 f("f", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
420 
421  f.SetParameter(0, h->GetMaximum());
422  f.SetParameter(1, h->GetMean());
423  f.SetParameter(2, h->GetRMS() * 0.25);
424  f.SetParameter(3, h->GetMinimum());
425 
426  h->Fit(&f, "Q", "same");
427 
428  bg = f.GetParameter(3) * h->GetNbinsX();
429  sg = h->GetSumOfWeights() - bg;
430  er = f.GetParError(3) * h->GetNbinsX();
431 
432  }
433 
434  int i = bin - 1;
435 
436  out[0][i] = sg;
437  out[1][i] = bg;
438  out[2][i] = er;
439 
440  delete h;
441 
442  }
443 
444  return out;
445 
446  }
447 
448 
449  /**
450  * Generate fitted L1 sample
451  */
452 
453  bg_type generate_fitted(int rb = 1) {
454 
455  bg_type out;
456 
457  int k = rnd->Integer(rt.size());
458 
459  int limit = rt[k].size() - (genRatio * TWindow_ms);
460 
461  int start = rnd->Integer(limit);
462 
463  h2d_t* sample = build_H2D(k, start);
464 
465  if (rb != 1) { sample->RebinX(rb); }
466 
467  bg_type fit = fit_H2D(sample);
468 
469  out.push_back(fit[0]);
470 
471  delete sample;
472 
473  if (rb == 1) {
474  out.push_back(build_NC(k, start));
475  } else {
476  out.push_back(rebin(build_NC(k, start), rb));
477  }
478 
479  out.push_back(fit[1]);
480  out.push_back(fit[2]);
481 
482  return out;
483 
484  }
485 
486  };
487 }
488 
489 #endif
JSUPERNOVA::JL0BackgroundSimulator::JL0BackgroundSimulator
JL0BackgroundSimulator(TFile *in, bool twoDim=false)
Default constructor.
Definition: JL0BackgroundSimulator.hh:83
JSUPERNOVA::JL0BackgroundSimulator::fit_H2D
static bg_type fit_H2D(h2d_t *in)
Fit 2D sample.
Definition: JL0BackgroundSimulator.hh:401
JTOOLS::n
const int n
Definition: JPolint.hh:628
JSUPERNOVA::JL0BackgroundSimulator
Class to emulate L0 background for an arbitrarily sized detector.
Definition: JL0BackgroundSimulator.hh:57
std::vector
Definition: JSTDTypes.hh:12
JSUPERNOVA::JL0BackgroundSimulator::genScale
double genScale
Definition: JL0BackgroundSimulator.hh:69
JTOOLS::j
int j
Definition: JPolint.hh:634
JSUPERNOVA::loadHistogram
vector< double > loadHistogram(TH1D *in)
Load histogram values to vector, each bin is converted to an element.
Definition: JL0BackgroundSimulator.hh:41
JSUPERNOVA::JL0BackgroundSimulator::~JL0BackgroundSimulator
~JL0BackgroundSimulator()
Destructor.
Definition: JL0BackgroundSimulator.hh:152
JSUPERNOVA::JL0BackgroundSimulator::genRatio
int genRatio
Definition: JL0BackgroundSimulator.hh:65
JSUPERNOVA::JL0BackgroundSimulator::build_NC
vector< double > build_NC(const int k, const int start)
Build NC sequence.
Definition: JL0BackgroundSimulator.hh:336
JSUPERNOVA::JL0BackgroundSimulator::nc
vector< vector< double > > nc
Definition: JL0BackgroundSimulator.hh:60
JSUPERNOVA::JL0BackgroundSimulator::rt2d
vector< h2d_t * > rt2d
Definition: JL0BackgroundSimulator.hh:62
JSUPERNOVA::JL0BackgroundSimulator::generate_fitted
bg_type generate_fitted(int rb=1)
Generate fitted L1 sample.
Definition: JL0BackgroundSimulator.hh:453
JSUPERNOVA::JL0BackgroundSimulator::generate_H2D
h2d_t * generate_H2D(int rb=1)
Generate 2D sample.
Definition: JL0BackgroundSimulator.hh:375
JSUPERNOVA::JL0BackgroundSimulator::rnd
TRandom * rnd
Definition: JL0BackgroundSimulator.hh:71
JSUPERNOVA::bg_type
vector< vector< double > > bg_type
Definition: JL0BackgroundSimulator.hh:32
JSUPERNOVA::JL0BackgroundSimulator::TWindow_ms
int TWindow_ms
Definition: JL0BackgroundSimulator.hh:64
JSUPERNOVA::JL0BackgroundSimulator::generate
bg_type generate()
Generate sample of L0 background L0 data are randomly sampled from a single L0 dataset.
Definition: JL0BackgroundSimulator.hh:204
JSUPERNOVA::JL0BackgroundSimulator::setSeed
void setSeed(const UInt_t uSeed=0)
Set TRandom seed.
Definition: JL0BackgroundSimulator.hh:170
JSUPERNOVA::JL0BackgroundSimulator::generate_poisson
bg_type generate_poisson(const int domRate_Hz=500)
Generate pure poissionian background.
Definition: JL0BackgroundSimulator.hh:233
std
Definition: jaanetDictionary.h:36
JSUPERNOVA::JL0BackgroundSimulator::configureRatio
void configureRatio(const int iSize, const int oSize)
Set rebinning ratio.
Definition: JL0BackgroundSimulator.hh:189
JSUPERNOVA::JL0BackgroundSimulator::configureTimeWindow
void configureTimeWindow(const int T_ms)
Configure the duration of an output sample.
Definition: JL0BackgroundSimulator.hh:161
JSUPERNOVA::JL0BackgroundSimulator::rt
vector< vector< double > > rt
Definition: JL0BackgroundSimulator.hh:59
JSUPERNOVA::JL0BackgroundSimulator::generate_shuffled
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 ...
Definition: JL0BackgroundSimulator.hh:258
JSUPERNOVA::JL0BackgroundSimulator::rebin
static vector< double > rebin(const vector< double > &in, const int rb)
Rebin vector.
Definition: JL0BackgroundSimulator.hh:358
JSUPERNOVA::JL0BackgroundSimulator::build_H2D
h2d_t * build_H2D(const int k, const int start)
Generate 2D histogram.
Definition: JL0BackgroundSimulator.hh:297
JSUPERNOVA::JL0BackgroundSimulator::outputSize
int outputSize
Definition: JL0BackgroundSimulator.hh:66
JSUPERNOVA::h2d_t
TH2F h2d_t
Definition: JRipple.hh:7
JSUPERNOVA
Definition: JSupernova.hh:24
JRipple.hh