---

DistRandomiser2D.cxx

Go to the documentation of this file.

#include "FastShowerUtils/Normalisers/DistRandomiser2D.h"

#include <cmath>
#include <string>
#include <vector>
#include <algorithm>
#include <iostream>
#include <fstream>
namespace FastShower{

  DistRandomiser2D::DistRandomiser2D(const std::string& f){

    // get data form file (histogram)
    std::vector<int> v;
    readData(f,v);
    
    // needed to initialise RandGeneral
    std::vector<int>::size_type nn =  v.size(); // # of xbins*ybins
    int n = (int) std::sqrt(nn); // # of bins
    // make sure vector has 50 rows, each with 50 elements
    assert(n == 50);

    // histograms
    double yDist[50], xDist[50], xyDist[50][50];

    // map vector into a matrix
    std::vector<int>::const_iterator itr = v.begin();
    std::vector<int>::const_iterator end = v.end();
    for (int j=0; j<n; ++j){
      for (int i=0; i<n; ++i){
        xyDist[i][j] = *itr;
        ++itr; 
        if (itr == end) itr=v.begin();
      }
    }
    // extract distributions
    for (int j=0; j<n; ++j) {
      xDist[j] = 0;
      for (int i=0; i<n; ++i){
        yDist[i] = 0;
        // projection of y onto x-axis
        xDist[j] += xyDist[i][j];
        // slice x-axis, project each x-slice onto y-axis
        yDist[i] = xyDist[i][j];
      }
      // setup distributions of xSlices
      m_xSlices.push_back(new RandGeneral(yDist, n, 0));
      //for (int nb=0; nb<n; ++nb) cout<<yDist[nb]<<" "; cout<<endl;
    }
    // setup distribution
    m_yProjection=new RandGeneral(xDist, n, 0);
    //for (int nb=0; nb<n; ++nb) cout<<xDist[nb]<<" "; cout<<endl;
  }
  //
  void
  DistRandomiser2D::readData(const std::string& f, std::vector<int>& v) const {
    ifstream data(f.c_str());
    istream_iterator<int> dataItr(data);
    istream_iterator<int> endOfData;
    copy(dataItr,endOfData,back_inserter(v));
  }
  //
  std::pair<double,double> DistRandomiser2D::sample() const {
    // sample Ecal fraction
    double eFrac = 2.0 * m_yProjection->fire();
    // map to bin number (gives Hcal distribution)
    int xSlice = (int) (eFrac/0.04); ++xSlice;
    // 
    std::vector<RandGeneral*>::const_iterator itr = m_xSlices.begin(); 
    for (int i=1; i<xSlice; ++i){++itr;}
    assert(itr<m_xSlices.end());
    // sample Hcal fraction
    double hFrac = 2.0 * (*itr)->fire();
    std::pair<double,double> ehFracs(eFrac,hFrac);
    return ehFracs;
  }
  //
  DistRandomiser2D::~DistRandomiser2D(){
  }
}//namespace

---