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 
68  double genScale;
69 
70  TRandom* rnd;
71 
72  public:
73 
74  /**
75  * Default constructor
76  *
77  * The input file contains a sequence of histograms with sampled L0 data;
78  * Each histogram is loaded as a vector, which will be considered as an individual L0 dataset.
79  *
80  * \param in input file from JRipple
81  * \param twoDim make 2D histogram with distribution of hit time differences vs time
82  */
83  JL0BackgroundSimulator(TFile* in, bool twoDim = false) {
84 
85  TWindow_ms = 100;
86  genRatio = 115;
87  genScale = 1;
88 
89  TString rt_tag = ".RT_DET_SUM";
90  TString nc_tag = ".NC_DET_SUM";
91  TString rt2d_tag = ".RT2D_DET";
92 
93  rnd = new TRandom();
94 
95  TIter iter(in->GetListOfKeys());
96 
97  for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
98 
99  const TString tag(key->GetName());
100 
101  // match
102  if (tag.EndsWith(rt_tag)) {
103 
104  TString rt_hist_tag = tag;
105  TString run_tag = TString(tag, tag.Length() - rt_tag.Length());
106 
107  TString nc_hist_tag = run_tag + nc_tag;
108 
109  TH1D* RT = (TH1D*) in->Get(rt_hist_tag);
110  TH1D* NC = (TH1D*) in->Get(nc_hist_tag);
111 
112  vector<double> rt_buf = loadHistogram(RT);
113  vector<double> nc_buf = loadHistogram(NC);
114 
115  int n = nc_buf.size();
116 
117  rt.push_back( vector<double>() );
118  nc.push_back( vector<double>() );
119 
120  // build parallel vectors
121  // 1 NC bin <=> 100 RT bins
122  for (int i = 0; i < n; i++) {
123  // empty slices are discarded
124  if (nc_buf[i] > 0) {
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] );
128  }
129  }
130  }
131 
132  if (twoDim) {
133  TString rt2d_hist_tag = run_tag + rt2d_tag;
134  h2d_t* RT2D = (h2d_t*) in->Get(rt2d_hist_tag);
135 
136  if (RT2D != NULL) {
137  RT2D->SetDirectory(0);
138  rt2d.push_back(RT2D);
139  }
140  }
141 
142  }
143 
144  }
145 
146  }
147 
148  /**
149  * Destructor
150  */
152  delete rnd;
153  }
154 
155  /**
156  * Configure the duration of an output sample
157  *
158  * \param T_ms sample duration in ms
159  */
160  void configureTimeWindow(const int T_ms) {
161  TWindow_ms = T_ms;
162  }
163 
164  /**
165  * Set TRandom seed
166  *
167  * \param uSeed seed
168  */
169  void setSeed(const UInt_t uSeed = 0) {
170  rnd->SetSeed(uSeed);
171  }
172 
173  /**
174  * Configure generation ratio
175  *
176  * \param inputNumberOfLines number of lines of input detector
177  * \param outputNumberOfLines number of lines of output detector
178  */
179  void configureRatio(const int inputNumberOfLines, const int outputNumberOfLines) {
180  outputSize = inputNumberOfLines;
181  // ratio is integer by design
182  genRatio = outputNumberOfLines / inputNumberOfLines;
183  // scaling factor to compensate the ratio remainder
184  genScale = outputNumberOfLines / (1.0 * inputNumberOfLines * genRatio);
185  }
186 
187 
188  /**
189  * Generate sample of L0 background
190  * L0 data are randomly sampled from a single L0 dataset.
191  *
192  * \return 2D vector with hit count and number of active channels as a function of time
193  */
195 
196  bg_type out(2, vector<double>(TWindow_ms));
197 
198  int k = rnd->Integer(rt.size());
199 
200  int limit = rt[k].size() - (genRatio * TWindow_ms);
201 
202  int start = rnd->Integer(limit);
203 
204  for (int i = 0; i < genRatio; i++) {
205 
206  for (int j = 0; j < TWindow_ms; j++) {
207 
208  int bin = start + (i * TWindow_ms) + j;
209 
210  out[0][j] += genScale * rt[k].at(bin );
211  out[1][j] += genScale * nc[k].at(bin / 100);
212 
213  }
214 
215  }
216 
217  return out;
218  }
219 
220  /**
221  * Generate pure poissionian background
222  */
223  bg_type generate_poisson(const int domRate_Hz = 500) {
224 
225  bg_type out(2, vector<double>(TWindow_ms));
226 
227  double mu = outputSize * 18 * domRate_Hz * 1.0e-3;
228 
229  for (unsigned i = 0; i < out[0].size(); i++) {
230 
231  double count = rnd->Poisson(mu);
232 
233  out[0][i] = count;
234  out[1][i] = outputSize * 18 * 31;
235 
236  }
237 
238  return out;
239 
240  }
241 
242 
243  /**
244  * Generate sample of L0 background
245  * The sampling of the L0 data is not sequential but random within the L0 dataset
246  * \return 2D vector with hit count and number of active channels as a function of time
247  */
248  bg_type generate_shuffled(bool randomizeRun = false) {
249 
250  bg_type out(2, vector<double>(TWindow_ms));
251 
252  int nRuns = rt.size();
253 
254  int k = rnd->Integer(nRuns);
255 
256  // a series n = 'genRatio' elements is sampled randomly within the dataset
257  // all the elements start with the same offset to the timeslice beginning
258 
259  int offset = rnd->Integer(100);
260 
261  for (int i = 0; i < genRatio; i++) {
262 
263  int start = 100 * rnd->Integer(nc[k].size() - ceil(TWindow_ms / 100.0)) + offset;
264 
265  for (int j = 0; j < TWindow_ms; j++) {
266 
267  int bin = start + j;
268 
269  out[0][j] += genScale * rt[k].at(bin );
270  out[1][j] += genScale * nc[k].at(bin / 100);
271 
272  }
273 
274  if (randomizeRun) {
275  k = rnd->Integer(nRuns);
276  }
277 
278  }
279 
280  return out;
281  }
282 
283 
284  /**
285  * Generate 2D histogram
286  */
287  h2d_t* build_H2D (const int k, const int start) {
288 
289  h2d_t* proto = rt2d[k];
290 
291  int ny = proto->GetNbinsY();
292  int ymin = proto->GetYaxis()->GetXmin();
293  int ymax = proto->GetYaxis()->GetXmax();
294 
295  h2d_t* out = new h2d_t("sample", NULL, TWindow_ms, 0, TWindow_ms, ny, ymin, ymax);
296 
297  for (int i = 0; i < genRatio; i++) {
298 
299  for (int j = 0; j < TWindow_ms; j++) {
300 
301  int xbin = 1 + start + (i * TWindow_ms) + j;
302 
303  for (int ybin = 1; ybin <= rt2d[k]->GetNbinsY(); ybin++) {
304 
305  double val = rt2d[k]->GetBinContent(xbin, ybin);
306 
307  double bcx = ((TAxis*)out->GetXaxis())->GetBinCenter(j + 1);
308  double bcy = ((TAxis*)out->GetYaxis())->GetBinCenter(ybin);
309 
310  out->Fill(bcx, bcy, val);
311 
312  }
313 
314  }
315 
316  }
317 
318  return out;
319 
320  }
321 
322 
323  /**
324  * Build NC sequence
325  */
326  vector<double> build_NC (const int k, const int start) {
327 
328  vector<double> out(TWindow_ms);
329 
330  for (int i = 0; i < genRatio; i++) {
331 
332  for (int j = 0; j < TWindow_ms; j++) {
333 
334  int bin = start + (i * TWindow_ms) + j;
335 
336  out[j] += genScale * nc[k].at(bin / 100);
337 
338  }
339 
340  }
341 
342  return out;
343  }
344 
345  /**
346  * Rebin vector
347  */
348  static vector<double> rebin(const vector<double> &in, const int rb) {
349 
350  const int n = in.size();
351 
352  vector<double> out(n/rb);
353 
354  for (int i = 0; i < n; i++) {
355  out[i/rb] += in[i];
356  }
357 
358  return out;
359 
360  }
361 
362  /**
363  * Generate 2D sample
364  */
365  h2d_t* generate_H2D(int rb = 1) {
366 
367  h2d_t* out = NULL;
368 
369  if (rt2d.size() > 0) {
370 
371  int k = rnd->Integer(rt2d.size());
372 
373  int limit = rt2d[k]->GetNbinsX() - (genRatio * TWindow_ms);
374 
375  int start = rnd->Integer(limit);
376 
377  out = build_H2D(k, start);
378 
379  }
380 
381  if (rb != 1) { out->RebinX(rb); }
382 
383  return out;
384 
385  }
386 
387 
388  /**
389  * Fit 2D sample
390  */
391  static bg_type fit_H2D(h2d_t* in) {
392 
393  int n = in->GetNbinsX();
394 
395  // 0 -> signal
396  // 1 -> fitted background
397  // 2 -> background fit error
398 
399  bg_type out(3, vector<double>(n));
400 
401  for (int bin = 1; bin <= n; bin++) {
402 
403  double sg = 0, bg = 0, er = 0;
404 
405  TH1D* h = in->ProjectionY("timeBin", bin, bin);
406 
407  if (h->GetEntries() > 0) {
408 
409  TF1 f("f", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
410 
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());
415 
416  h->Fit(&f, "Q", "same");
417 
418  bg = f.GetParameter(3) * h->GetNbinsX();
419  sg = h->GetSumOfWeights() - bg;
420  er = f.GetParError(3) * h->GetNbinsX();
421 
422  }
423 
424  int i = bin - 1;
425 
426  out[0][i] = sg;
427  out[1][i] = bg;
428  out[2][i] = er;
429 
430  delete h;
431 
432  }
433 
434  return out;
435 
436  }
437 
438 
439  /**
440  * Generate fitted L1 sample
441  */
442 
443  bg_type generate_fitted(int rb = 1) {
444 
445  bg_type out;
446 
447  int k = rnd->Integer(rt.size());
448 
449  int limit = rt[k].size() - (genRatio * TWindow_ms);
450 
451  int start = rnd->Integer(limit);
452 
453  h2d_t* sample = build_H2D(k, start);
454 
455  if (rb != 1) { sample->RebinX(rb); }
456 
457  bg_type fit = fit_H2D(sample);
458 
459  out.push_back(fit[0]);
460 
461  delete sample;
462 
463  if (rb == 1) {
464  out.push_back(build_NC(k, start));
465  } else {
466  out.push_back(rebin(build_NC(k, start), rb));
467  }
468 
469  out.push_back(fit[1]);
470  out.push_back(fit[2]);
471 
472  return out;
473 
474  }
475 
476  };
477 }
478 
479 #endif
*fatal Wrong number of arguments esac JCookie sh JRuns D $DETECTOR d sort n k
Definition: JRunrange.sh:16
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.
esac print_variable DETECTOR INPUT_FILE OUTPUT_FILE CDF for TYPE in
Definition: JSirene.sh:45
vector< vector< double > > bg_type
h2d_t * build_H2D(const int k, const int start)
Generate 2D histogram.
then JPizza f
Definition: JPizza.sh:46
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 ...
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.
void configureRatio(const int inputNumberOfLines, const int outputNumberOfLines)
Configure generation ratio.
static vector< double > rebin(const vector< double > &in, const int rb)
Rebin vector.
std::vector< int > count
Definition: JAlgorithm.hh:184
alias put_queue eval echo n
Definition: qlib.csh:19
void configureTimeWindow(const int T_ms)
Configure the duration of an output sample.
int j
Definition: JPolint.hh:634
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.