JLightCurveBackgroundGenerator.hh

Go to the documentation of this file.

#ifndef __JSUPERNOVA_JLIGHTCURVEBACKGROUNDGENERATOR__
#define __JSUPERNOVA_JLIGHTCURVEBACKGROUNDGENERATOR__
 
#include <vector>
 
#include "TFile.h"
#include "TKey.h"
 
#include "TH1D.h"
#include "TH2F.h"
#include "TF1.h"
 
#include "TRandom3.h"
 
#include "JROOT/JMinimizer.hh"
 
#include "JRipple.hh"
 
/**
 * \file
 * 
 * This file provides a class to emulate a L0 (L1) background on the ms time scale for an arbitrarily sized detector.
 * The emulation is based on the sampling of real data (in principle from a small-sized detector) through the JRipple application.
 * Samples of given duration are taken in a contiguous sequence and summed bin-by-bin as they were coming from different DUs,
 * in this way time-correlations are used to emulate space-correlations which are not possible to probe due to the detector limited size.
 *
 * \author mlincett
 */
 
namespace JSUPERNOVA {
 
  using namespace std;
 
  typedef vector<vector <double> > bg_type; 
 
  /**
   * Load histogram values to vector, each bin is converted to an element.
   * No other information is preserved (binning, axis semantics, under/overflows).
   *
   * \param  in           input histogram
   * \return              vector with one value per bin              
   */
  vector<double> loadHistogram(TH1D* in) {
    int n = in->GetNbinsX();
    
    vector<double> out(n);
    
    for (int i = 0; i < n; i++) {
      out[i] = in->GetBinContent(i + 1);
    }
 
    return out;
  }
 
  /**
   * Class to emulate L0 background for an arbitrarily sized detector.
   */
 
  class JLightCurveBackgroundGenerator {
  private:
    vector<vector<double> > rt;
    vector<vector<double> > nc;
 
    vector<h2d_t*> rt2d;
 
    int    TWindow_ms;
    int    genRatio;  
    int    outputSize;
 
    double genScale;  
 
    TRandom* rnd;
 
  public:
 
    /**
     * Default constructor
     *
     * The input file contains a sequence of histograms with sampled L0 data;
     * Each histogram is loaded as a vector, which will be considered as an individual L0 dataset.
     *
     * \param  in           input file from JRipple
     * \param twoDim        make 2D histogram with distribution of hit time differences vs time
     */
    JLightCurveBackgroundGenerator(TFile* in, bool twoDim = false) {
 
      TWindow_ms = 100;
      genRatio   = 115;
      genScale   = 1;
 
      TString rt_tag = ".RT_DET_SUM";
      TString nc_tag = ".NC_DET_SUM";
      TString rt2d_tag = ".RT2D_DET";
 
      rnd = new TRandom();
 
      TIter iter(in->GetListOfKeys());
 
      for (TKey* key; (key = (TKey*) iter.Next()) != NULL; ) {
 
        const TString tag(key->GetName());
 
        // match 
        if (tag.EndsWith(rt_tag)) {
 
          TString rt_hist_tag = tag;
          TString run_tag = TString(tag, tag.Length() - rt_tag.Length());
          
          TString nc_hist_tag = run_tag + nc_tag;
 
          TH1D* RT = (TH1D*) in->Get(rt_hist_tag);
          TH1D* NC = (TH1D*) in->Get(nc_hist_tag);
 
          vector<double> rt_buf = loadHistogram(RT);
          vector<double> nc_buf = loadHistogram(NC);
 
          int n = nc_buf.size();
 
          rt.push_back( vector<double>() );
          nc.push_back( vector<double>() );
 
          // build parallel vectors 
          // 1 NC bin <=> 100 RT bins
          for (int i = 0; i < n; i++) {
            // empty slices are discarded
            if (nc_buf[i] > 0) {
              nc.back().push_back( nc_buf[i] );
              for (int j = 100 * i; j < 100 * (i + 1); j++) {
                rt.back().push_back( rt_buf[j] );
              }
            }
          }
        
          if (twoDim) {
            TString rt2d_hist_tag = run_tag + rt2d_tag;
            h2d_t* RT2D = (h2d_t*) in->Get(rt2d_hist_tag);
 
            if (RT2D != NULL) {
              RT2D->SetDirectory(0);
              rt2d.push_back(RT2D);
            }
          }
 
        }
 
      }
 
    }
 
    /**
     * Destructor
     */
    ~JLightCurveBackgroundGenerator() {
      delete rnd;
    }
    
    /**
     * Configure the duration of an output sample
     *
     * \param  T_ms sample duration in ms
     */
    void configureTimeWindow(const int T_ms) {
      TWindow_ms = T_ms;
    }
 
    /**
     * Set TRandom seed
     *
     * \param  uSeed seed
     */
    void setSeed(const UInt_t uSeed = 0) {
      rnd->SetSeed(uSeed);
    } 
    
    /**
     * Configure generation ratio
     *
     * \param  inputNumberOfLines  number of lines of input detector
     * \param  outputNumberOfLines  number of lines of output detector
     */
    void configureRatio(const int inputNumberOfLines, const int outputNumberOfLines) {
      outputSize   = inputNumberOfLines;
      // ratio is integer by design
      genRatio     = outputNumberOfLines / inputNumberOfLines;
      // scaling factor to compensate the ratio remainder
      genScale     = outputNumberOfLines / (1.0 * inputNumberOfLines * genRatio);
    }
    
 
    /**
     * Generate sample of L0 background
     * L0 data are randomly sampled from a single L0 dataset.
     *
     * \return 2D vector with hit count and number of active channels as a function of time
     */
    bg_type generate() {
 
      bg_type out(2, vector<double>(TWindow_ms));
 
      int k = rnd->Integer(rt.size());
 
      int limit = rt[k].size() - (genRatio * TWindow_ms);
 
      int start = rnd->Integer(limit);
 
      for (int i = 0; i < genRatio; i++) {
 
        for (int j = 0; j < TWindow_ms; j++) {
 
          int bin = start + (i * TWindow_ms) +  j;
 
          out[0][j] += genScale * rt[k].at(bin      );
          out[1][j] += genScale * nc[k].at(bin / 100);
 
        }
 
      }
 
      return out;
    }
 
    /** 
     * Generate pure poissionian background
     */
    bg_type generate_poisson(const int domRate_Hz = 500) {
 
      bg_type out(2, vector<double>(TWindow_ms));
 
      double mu = outputSize * 18 * domRate_Hz * 1.0e-3; 
 
      for (unsigned i = 0; i < out[0].size(); i++) {
 
        double count  = rnd->Poisson(mu);
        
        out[0][i] = count;
        out[1][i] = outputSize * 18 * 31;
        
      } 
 
      return out;
 
    }
 
    
    /**
     * Generate sample of L0 background
     * The sampling of the L0 data is not sequential but random within the L0 dataset
     * \return 2D vector with hit count and number of active channels as a function of time
     */
    bg_type generate_shuffled(bool randomizeRun = false) {
 
      bg_type out(2, vector<double>(TWindow_ms));
 
      int nRuns = rt.size();
 
      int k = rnd->Integer(nRuns);
 
      // a series n = 'genRatio' elements is sampled randomly within the dataset
      // all the elements start with the same offset to the timeslice beginning
 
      int offset = rnd->Integer(100);
 
      for (int i = 0; i < genRatio; i++) {
 
        int start = 100 * rnd->Integer(nc[k].size() - ceil(TWindow_ms / 100.0)) + offset;
 
        for (int j = 0; j < TWindow_ms; j++) {
 
          int bin = start + j;
 
          out[0][j] += genScale * rt[k].at(bin      );
          out[1][j] += genScale * nc[k].at(bin / 100);
                  
        }
 
        if (randomizeRun) {
          k = rnd->Integer(nRuns);
        }
 
      }
 
      return out;
    }    
    
    
    /**
     * Generate 2D histogram
     */
    h2d_t* build_H2D (const int k, const int start) {
 
      h2d_t* proto = rt2d[k];
 
      int ny   = proto->GetNbinsY();
      int ymin = proto->GetYaxis()->GetXmin(); 
      int ymax = proto->GetYaxis()->GetXmax();
 
      h2d_t* out = new h2d_t("sample", NULL, TWindow_ms, 0, TWindow_ms, ny, ymin, ymax);        
      
      for (int i = 0; i < genRatio; i++) {
 
        for (int j = 0; j < TWindow_ms; j++) {
 
          int xbin = 1 + start + (i * TWindow_ms) +  j;
 
          for (int ybin = 1; ybin <= rt2d[k]->GetNbinsY(); ybin++) {
              
            double val = rt2d[k]->GetBinContent(xbin, ybin);    
 
            double bcx = ((TAxis*)out->GetXaxis())->GetBinCenter(j + 1);
            double bcy = ((TAxis*)out->GetYaxis())->GetBinCenter(ybin);
              
            out->Fill(bcx, bcy, val);
 
          }
 
        }
 
      }
 
      return out;
 
    }
 
 
    /**
     * Build NC sequence
     */
    vector<double> build_NC (const int k, const int start) {
      
      vector<double> out(TWindow_ms);
 
      for (int i = 0; i < genRatio; i++) {
 
        for (int j = 0; j < TWindow_ms; j++) {
 
          int bin = start + (i * TWindow_ms) +  j;
 
          out[j] += genScale * nc[k].at(bin / 100);
 
        }
 
      }
 
      return out;
    }
 
    /** 
     * Rebin vector
     */
    static vector<double> rebin(const vector<double> &in, const int rb) {
 
      const int n = in.size();
 
      vector<double> out(n/rb);
     
      for (int i = 0; i < n; i++) {
        out[i/rb] += in[i];
      }
 
      return out;
 
    }
 
    /**
     * Generate 2D sample
     */
    h2d_t* generate_H2D(int rb = 1) {
 
      h2d_t* out = NULL;
      
      if (rt2d.size() > 0) {
 
        int k = rnd->Integer(rt2d.size());
 
        int limit = rt2d[k]->GetNbinsX() - (genRatio * TWindow_ms);
 
        int start = rnd->Integer(limit);
 
        out = build_H2D(k, start);
 
      }
 
      if (rb != 1) { out->RebinX(rb); }
 
      return out;
    
    }
 
    
    /**
     * Fit 2D sample
     */
    static bg_type fit_H2D(h2d_t* in) {
 
      int n = in->GetNbinsX();
 
      // 0 -> signal
      // 1 -> fitted background
      // 2 -> background fit error
 
      bg_type out(3, vector<double>(n));
 
      for (int bin = 1; bin <= n; bin++) {
 
        double sg = 0, bg = 0, er = 0;
 
        TH1D* h = in->ProjectionY("timeBin", bin, bin);
 
        if (h->GetEntries() > 0) {
          
          TF1 f("f", "[0]*exp(-0.5*(x-[1])*(x-[1])/([2]*[2]))/(TMath::Sqrt(2*TMath::Pi())*[2]) + [3]");
 
          f.SetParameter(0, h->GetMaximum());
          f.SetParameter(1, h->GetMean());
          f.SetParameter(2, h->GetRMS() * 0.25);
          f.SetParameter(3, h->GetMinimum());
          
          h->Fit(&f, "Q", "same");
 
          bg = f.GetParameter(3) * h->GetNbinsX();
          sg = h->GetSumOfWeights() - bg;
          er = f.GetParError(3) * h->GetNbinsX();
 
        }
 
        int i = bin - 1;
      
        out[0][i] = sg;
        out[1][i] = bg;
        out[2][i] = er;
 
        delete h;
 
      }
      
      return out;
 
    }
 
 
    /**
     * Generate fitted L1 sample 
     */
 
    bg_type generate_fitted(int rb = 1) {
 
      bg_type out;
 
      int k = rnd->Integer(rt.size());
 
      int limit = rt[k].size() - (genRatio * TWindow_ms);
 
      int start = rnd->Integer(limit);
 
      h2d_t* sample = build_H2D(k, start);
      
      if (rb != 1) { sample->RebinX(rb); }
 
      bg_type fit = fit_H2D(sample);
      
      out.push_back(fit[0]);
 
      delete sample;
      
      if (rb == 1) {
        out.push_back(build_NC(k, start));
      } else {
        out.push_back(rebin(build_NC(k, start), rb));
      }
 
      out.push_back(fit[1]);
      out.push_back(fit[2]);
 
      return out;      
    
    }   
 
  };
}
 
#endif