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