JL0BackgroundSimulator.hh

Go to the documentation of this file.

#ifndef __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__
#define __JSUPERNOVA_JL0BACKGROUNDSIMULATOR__

#include <vector>

#include "TFile.h"
#include "TKey.h"

#include "TH1D.h"
#include "TH2F.h"
#include "TF1.h"

#include "TRandom3.h"

#include "JRipple.hh"

/**
 * \file
 * 
 * This file provides a class to emulate an L0 background on the ms time scale for an arbitrarily sized detector.
 * The emulation is based on the sampling of real L0 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 JL0BackgroundSimulator {
  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
     */
    JL0BackgroundSimulator(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
     */
    ~JL0BackgroundSimulator() {
      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