RawScanResultStreamer_v3.cpp

#include "RawScanResultStreamer_v3.h"
#include "Sct/SctNames.h"
#include "Sct/SctParameters.h"
#include "Sct/LogicErrors.h"
#include "Sct/VersionNotSupportedException.h"
#include "Sct/UnsupportedOperationError.h"
#include "../RawScanResult.h"
#include "../ConfigurationVariable.h"
#include "ScanResultWriter/dataTypes.h"

#include "TH2.h"

#include <iostream>
#include <sstream>
#include <boost/lexical_cast.hpp>

using namespace std;
using namespace Sct;
using namespace boost;

namespace SctData {
namespace IO {

/* READ NOTES ON STREAMERS AND VERSIONS BEFORE EDITING THIS FILE! */
unsigned RawScanResultStreamer_v3::s_version=3;

RawScanResultStreamer_v3::RawScanResultStreamer_v3() throw() {}

bool RawScanResultStreamer_v3::inMap = IOManager::addToMap("SctData::RawScanResult",  auto_ptr<Streamer>(new RawScanResultStreamer_v3()));

shared_ptr<Streamable> RawScanResultStreamer_v3::read(IStream& in, const IOManager& manager) const {
    shared_ptr<Streamable> ad (&helper.create());
    read(in, *ad, manager);
    return ad;
}

void RawScanResultStreamer_v3::write(OStream& out, const Streamable& ob, const IOManager& manager) const {
    manager.writeImpl(out, ob, "SctData::ScanResult");

    //Downcast should always work!
    const RawScanResult& raw = dynamic_cast<const RawScanResult&>(ob);
    helper.set(raw);

    //Now publish data - we publish using the ROOTHIST format.
    TH2D* link0 = helper.getScanData(0);
    TH2D* link1 = helper.getScanData(1);
    out << (link0->GetNbinsX()+2)*(link0->GetNbinsY()+2);   //Note, size is number in 1 link ie no bins * no points
    out << SR_DT_ROOTHIST << SR_WD_64;
    //cout << "Raw.getDataType() is " << raw.getDataType() << endl;
    // ROOTDATA-specific stuff:
    out << raw.getDataType();
    out.put ( link0->GetArray(), (link0->GetNbinsX()+2)*(link0->GetNbinsY()+2) );
    out.put ( link1->GetArray(), (link1->GetNbinsX()+2)*(link1->GetNbinsY()+2) );

    unsigned nOcc = raw.nOccupancyPerEvent();
    out << (unsigned) nOcc;
    for (unsigned iOcc=0; iOcc<nOcc; ++iOcc){
      shared_ptr<TH2D> occ = helper.getOccupancyPerEvent(iOcc);
      if (occ.get()==0) {
    ostringstream oss;
    oss << "No occupancy histogram for chip" << iOcc << " of " << nOcc;
    throw LogicError(oss.str(), __FILE__, __LINE__ );
      } 
      out.put( occ->GetArray(), (occ->GetNbinsX()+2)*(occ->GetNbinsY()+2) );
    }
}

void RawScanResultStreamer_v3::read(IStream& in, Streamable& ob, const IOManager& manager) const {
    manager.readImpl(in, ob, "SctData::ScanResult");

    RawScanResult& raw = dynamic_cast<RawScanResult&>(ob);
    helper.set(raw);

    //cout <<"Parent finished" << endl;
    //cout << raw.getHeader().getRunNumber() << "  " << raw.getHeader().getScanNumber()
    //     << "  " << raw.getHeader().getModuleName() << endl;

    //Now refresh and fill the data
    readData(in, raw);
    
    //cout <<"Done" << endl;
}


void RawScanResultStreamer_v3::setHistSize(RawScanResult& raw, unsigned nbinsx, short unsigned ilink) const{
  // check its not already the right size!
  if ( helper.getScanData(ilink) && (unsigned) helper.getScanData(ilink)->GetNbinsX()==nbinsx ) return;

  //Create hists
  double* bins = raw.getPoints().getEdgesAscending();
  ostringstream name; name << (string)raw.getUniqueID() << "_link" << ilink;
  auto_ptr<TH2D> link_hist (new TH2D(name.str().c_str(), name.str().c_str(), nbinsx, -0.5, nbinsx-0.5,
                     raw.getPoints().getNPoints(), bins));
  
  delete [] bins;
  
  link_hist->SetEntries(1);
  helper.setScanData(ilink, link_hist);
}

void RawScanResultStreamer_v3::readData(IStream& in, RawScanResult& raw) const {
    unsigned int size;
    unsigned short type;
    unsigned short width;
    in >> size >> type >> width;
    helper.setDataType(type);

    switch (type) {
    case SR_DT_SLICE:
        readSliceData(size, width, in, raw);
        return;
    case SR_DT_ROOTHIST:
        readRootData(size, width, in, raw);
        return;
    case SR_DT_RAWHIST:
        readRawData(size, width, in, raw);
        return;
    default:
    ostringstream oss;
    oss << "RawScanResultIS::refreshData Unknown data format: " << type << " size was: " << size;
        throw VersionNotSupportedException(oss.str(), __FILE__, __LINE__); 
    }
}

void RawScanResultStreamer_v3::readRootData(unsigned int isize, unsigned short width, IStream& in, RawScanResult& raw) const {

  unsigned short original_type;
  in >> original_type;
  helper.setDataType(original_type);

    for (unsigned ilink=0; ilink<2; ++ilink){
      double* data;
      unsigned int size = 0;
      in.get(&data, size);    
      
      unsigned size_y= raw.getPoints().getNPoints() + 2 ;
      if (size % size_y !=0 ) {
    ostringstream oss;
    oss << "RawScanResultIS::readRootData.  Link " << ilink << " size = " << size << ", not divisible by " << size_y;
    throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
      }
      unsigned size_x=size/size_y;
      if (size_x<=1){
    ostringstream oss;
    oss << "X-size of root data must be 2 or greater! Got "<< size_x;
    throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
      }
      setHistSize(raw, size_x-2, ilink);
      for (unsigned int i=0; i<size; ++i) {
        helper.getScanData(ilink)->SetBinContent(i, data[i]);
      }
      delete [] data;
    }

    unsigned nOcc = 0;  in >> nOcc;
    helper.setNumberOccupancyHists(nOcc);
    for (unsigned iOcc=0; iOcc<nOcc; ++iOcc){
      createOccupancyHistogram(helper, raw, iOcc);
      shared_ptr<TH2D> occ = helper.getOccupancyPerEvent(iOcc);
      if (!occ.get()) throw LogicError("Failed to create occupancy histogram", __FILE__, __LINE__);
      double* occ_data;
      unsigned int size = 0;
      in.get(&occ_data, size);
      
      const unsigned expectedSize=(occ->GetNbinsX()+2)*(occ->GetNbinsY()+2);

      if ( size!= expectedSize ){
    ostringstream oss;
    oss << "Occupancy Histogram data size was " << size
        << " expected " << expectedSize << " (Chip " << iOcc<<")";
    throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
      }
      for (unsigned ibin=0; ibin<size; ++ibin){
    occ->SetBinContent(ibin, occ_data[ibin]);
      }
      delete [] occ_data;
    }
}

void RawScanResultStreamer_v3::readRawData(unsigned int isize, unsigned short width, IStream& in, RawScanResult& raw) const {
  int npoints = raw.getPoints().getNPoints();
  unsigned wid_mult = width==SR_WD_16 ? 1 : width==SR_WD_32 ? 2 : width==SR_WD_64 ? 4 : 0;
  unsigned quotient = npoints * wid_mult;
  if (quotient==0 || isize % quotient != 0) {
    ostringstream oss; oss << "Data size  = " << isize << " but trying to divide by " << quotient;
    throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
  } 
  unsigned x_size = isize / quotient;

  for (unsigned ilink=0; ilink<2; ++ilink){
    setHistSize(raw, x_size, ilink);
  }
  
  TH2D* link0 = helper.getScanData(0);
  TH2D* link1 = helper.getScanData(1);

  unsigned size=0;
  bool ascending=raw.getPoints().ascending();
  switch (width) {
  case SR_WD_16:{
    UINT16* data = 0;
    in.get(&data, size);
    if (size!=isize){
      ostringstream oss; oss << "Size (" << size << ") dosent match header size(" << isize << ")";
      throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
    }
    for (int j=0; j<npoints; ++j) {
      int ipt=ascending ? j : npoints-j-1;
      for (unsigned int i=0; i<x_size; ++i) {
    //+1 because bin 0 is underflow
    link0->SetBinContent(link0->GetBin(i+1,j+1), data[ipt*x_size*2 + i]);
    link1->SetBinContent(link1->GetBin(i+1,j+1), data[ipt*x_size*2 + i+x_size]);
      }
    }
    delete[] data; return;
  }
  case SR_WD_32: {
    UINT32* data = 0;
    in.get(&data, size);
    if (size!=isize){
      ostringstream oss; oss << "Size (" << size << ") dosent match header size(" << isize << ")";
      throw StreamCorruptedException(oss.str(), __FILE__, __LINE__);
    }
    for (int j=0; j<npoints; ++j) {
      int ipt=ascending ? j : npoints-j-1;
      for (unsigned int i=0; i<x_size; ++i) {
    //+1 because bin 0 is underflow
    link0->SetBinContent(link0->GetBin(i+1,j+1), data[ipt*x_size*2 + i]);
    link1->SetBinContent(link1->GetBin(i+1,j+1), data[ipt*x_size*2 + i+x_size]);
      }
    }
    delete[] data; return;
  }
  default : {
    ostringstream oss; oss << "Dont know how to interpret data of width " << width;
    throw UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
  }
  }
}

void RawScanResultStreamer_v3::readSliceData(unsigned int isize, unsigned short width, IStream& in, RawScanResult& raw) const {
    int npoints = raw.getPoints().getNPoints();
    unsigned expectedSize=2048u * npoints;

    //cout << width << endl;
    if (isize != expectedSize) 
    throw StreamCorruptedException("RawScanResultIS::readSliceData. Expected and header sizes not equal.  Expected: " + lexical_cast<string>(expectedSize) + 
                       " but header: " + lexical_cast<string>(isize), __FILE__, __LINE__);
    
    setHistSize(raw, nChannelLink, 0);
    setHistSize(raw, nChannelLink, 1);
    
    TH2D* link0 = helper.getScanData(0);
    TH2D* link1 = helper.getScanData(1);
    unsigned int size = 0;
    bool ascending=raw.getPoints().ascending();

    UINT32* data = 0; // standard data type to convert to...

    // copy or reinterpret data as correct size:
    switch (width) {
      //------------------------------------------------------------------------
    case SR_WD_16:{
      UINT16* data16;
      in.get(&data16, size);
      data=new UINT32[size];
      for (unsigned i=0; i<size; ++i){
        data[i]=data16[i];  // copy
      }
      delete[] data16;
      break;
    }
    case SR_WD_32: {
      in.get(&data, size);
      break;
    }
    default : {
      ostringstream oss; oss << "Dont know how to interpret data of width " << width;
      throw UnsupportedOperationError(oss.str(), __FILE__, __LINE__);
    }
    } // end of switch
    
    if (size != expectedSize){
      delete[] data;
      ostringstream msg;
      msg << "RawScanResultIS::readSliceData. Expected and actual sizes not equal.  Expected: "
      << expectedSize << " but actual: "
      << size << "; (DataWidth="<<width<<")";
      throw StreamCorruptedException(msg.str(), __FILE__, __LINE__);
    }
    
    for (int j=0; j<npoints; ++j) {
      int ipt=ascending ? j : npoints-j-1;
      for (unsigned int i=0; i<nChannelLink; ++i) {
    //+1 because 5in 0 is underflow
    link0->SetBinContent(link0->GetBin(i+1,j+1), data[ipt*2048 + i]);
    link1->SetBinContent(link1->GetBin(i+1,j+1), data[ipt*2048 + i+1024]);
      }
    }
    // Also get Occupancy histogram information from the 
    // memory location where the absent chips "6,7", "14,15" would be
    try{
      helper.setNumberOccupancyHists(nChipModule);
      for (unsigned ichip=0; ichip<nChipModule; ++ichip){
    createOccupancyHistogram(helper, raw, ichip);
    shared_ptr<TH2D> occ = helper.getOccupancyPerEvent(ichip);
    if (occ.get()==0) {
      ostringstream oss;
      oss << "Occupancy histogram does not exist, module"
          << raw.getHeader().getModuleName() << " chip " << ichip << endl;
      throw LogicError(oss.str(), __FILE__, __LINE__);
    } 
    // loop over scan points
    for (int j=0; j<npoints; ++j) {
      int ipt=ascending ? j : npoints-j-1;
          // calculate empty events from number of events - number of events with any data
      // loop over histogram bins
          double empty_events= raw.getPoints().getNEvents(ipt);
      for (int ibin=1; ibin<occ->GetNbinsX(); ++ibin){
        unsigned word = getOccWordPosition(ichip, ibin-1, ipt);
        occ->SetBinContent(ibin+1, j+1, data[word]);
            empty_events-=data[word];
      }
      occ->SetBinContent(1, j+1, empty_events);
    }
    occ->SetEntries(1);
      }
    }catch(Sct::LogicError& e){
      e.sendToMrs();
    }
    delete[] data;
  }

  void RawScanResultStreamer_v3::createOccupancyHistogram(const RawScanResultIOHelper& helper, const RawScanResult& raw, unsigned ichip){

    // name
    std::ostringstream oss;
    oss << raw.getHeader().getModuleName() << "_OccEvent_Chip_" << ichip;
    const char* name=oss.str().c_str();

    // xbins
    const unsigned nEdgesX=1+128/4 + 1;
    Double_t edgesX[nEdgesX];
    edgesX[0] = -0.5;
    edgesX[1] =  0.5;
    for(int i=2; i<nEdgesX; i++) {
      edgesX[i] = (i-1)*4+0.5;
    }

    // ybins
    const SctData::ScanPoints& points = raw.getPoints();
    const double* edgesY = points.getEdgesAscending();

    // create
    auto_ptr<TH2D> occ ( new TH2D(name, name, nEdgesX-1, edgesX, points.getNPoints(), edgesY) );
    
    // axis titles
    occ->SetXTitle("Channel count");
    occ->SetYTitle(raw.getHeader().getVariable().getVariableName().c_str());
    
    // clean up
    delete[] edgesY;

    // check
    if (!occ.get()) {
      throw LogicError("Could not make occupancy histogram", __FILE__, __LINE__);
    }

    // set
    helper.setOccupancyPerEvent(ichip, occ);
  }

  unsigned RawScanResultStreamer_v3::getOccWordPosition(unsigned ichip, unsigned ibin, unsigned ipt){
    return 128*6                                     // skip 6 chips worth of data
      + (ichip/6)*1024                               // + offset for link
      + (ichip%6)*32                                 // + start of chip 
      + ibin                                         // + bin
      + 2048 * ipt;                                  // + offset for data point
  }

}
}

---