Jpp
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
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
JL0BackgroundSimulator(TFile *in, bool twoDim=false)
Default constructor.
TH2F h2d_t
Definition: JRipple.hh:7
bg_type generate_fitted(int rb=1)
Generate fitted L1 sample.
vector< vector< double > > bg_type
h2d_t * build_H2D(const int k, const int start)
Generate 2D histogram.
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 ...
void configureRatio(const int iSize, const int oSize)
Set rebinning ratio.
static bg_type fit_H2D(h2d_t *in)
Fit 2D sample.
Class to emulate L0 background for an arbitrarily sized detector.
void setSeed(const UInt_t uSeed=0)
Set TRandom seed.
vector< double > build_NC(const int k, const int start)
Build NC sequence.
static vector< double > rebin(const vector< double > &in, const int rb)
Rebin vector.
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.