// ResponseMatrix object code. Definitions in ResponseMatrix.h

#ifndef HAVE_rmf
#include "rmf.h"
#endif

#ifndef HAVE_arf
#include "arf.h"
#endif

#ifndef HAVE_pha
#include "pha.h"
#endif

#ifndef HAVE_grouping
#include "grouping.h"
#endif

#ifndef HAVE_SPio
#include "SPio.h"
#endif

#ifndef HAVE_SPutils
#include "SPutils.h"
#endif

// Class rmf

// default constructor

rmf::rmf()
  : m_FirstChannel(0),
    m_NumberGroups(),
    m_FirstChannelGroup(),
    m_NumberChannelsGroup(),
    m_OrderGroup(),
    m_LowEnergy(),
    m_HighEnergy(),
    m_Matrix(),
    m_ChannelLowEnergy(),
    m_ChannelHighEnergy(),
    m_AreaScaling(1.0),
    m_ResponseThreshold(0.0),
    m_EnergyUnits("keV"),
    m_RMFUnits(" "),
    m_ChannelType("PI"),
    m_Telescope(" "),
    m_Instrument(" "),
    m_Detector(" "),
    m_Filter(" "),
    m_RMFType(" "),
    m_RMFExtensionName("MATRIX"),
    m_EBDExtensionName("EBOUNDS")
{
}


// Constructor from file

rmf::rmf(const string filename, const Integer RMFnumber)
{
  Integer status = read(filename, RMFnumber);
  if ( status != 0 ) clear();
}

// constructor from input arrays

rmf::rmf(const vector<Integer>& inNumberGroups,
	 const vector<vector<Integer> >& inFirstChannelGroup,
	 const vector<vector<Integer> >& inNumberChannelsGroup,
	 const vector<vector<Integer> >& inOrderGroup,
	 const vector<Real>& inLowEnergy, const vector<Real>& inHighEnergy,
	 const vector<vector<Real> >& inMatrix, 
	 const vector<Real>& inChannelLowEnergy,
	 const vector<Real>& inChannelHighEnergy,
	 const Integer inFirstChannel, const Real inAreaScaling,
	 const Real inResponseThreshold, const map<string,string>& inKeys)
  : m_FirstChannel(0),
    m_NumberGroups(),
    m_FirstChannelGroup(),
    m_NumberChannelsGroup(),
    m_OrderGroup(),
    m_LowEnergy(),
    m_HighEnergy(),
    m_Matrix(),
    m_ChannelLowEnergy(),
    m_ChannelHighEnergy(),
    m_AreaScaling(1.0),
    m_ResponseThreshold(0.0),
    m_EnergyUnits("keV"),
    m_RMFUnits(" "),
    m_ChannelType("PI"),
    m_Telescope(" "),
    m_Instrument(" "),
    m_Detector(" "),
    m_Filter(" "),
    m_RMFType(" "),
    m_RMFExtensionName("MATRIX"),
    m_EBDExtensionName("EBOUNDS")
{
  Integer status = load(inNumberGroups, inFirstChannelGroup, inNumberChannelsGroup,
			inOrderGroup, inLowEnergy, inHighEnergy, inMatrix,
			inChannelLowEnergy, inChannelHighEnergy, inFirstChannel,
			inAreaScaling, inResponseThreshold, inKeys);
  if ( status != OK ) clear();
}


// reading Matrix and Channel bounds extensions from RMF file.

Integer rmf::read(const string filename, const Integer RMFnumber)
{
  Integer Status(OK);
  
  // read the MATRIX extension

  Status = this->readMatrix(filename, RMFnumber);
  if ( Status != OK ) return(Status);

  // try to read the EBOUNDS extension. If this fails try to read with RMFnumber=1
  // under assumption that there are multiple MATRIX extensions and only one 
  // EBOUNDS extension. Read the EBOUNDS extension into a temporary rmf object so
  // that we can test for consistency between the MATRIX and the EBOUNDS entries

  rmf eBounds;

  Status = eBounds.readChannelBounds(filename, RMFnumber);
  if ( Status != OK ) {
    Status = eBounds.readChannelBounds(filename, 1);
  }

  // check for consistency between this and eBounds.

  if ( m_ChannelType != eBounds.getChannelType() ) {
    string msg = "CHANTYPE keyword differs between MATRIX and EBOUNDS extensions.";
    SPreportError(InconsistentKeywordValues, msg);
  }
  if ( m_Telescope != eBounds.getTelescope() ) {
    string msg = "TELESCOP keyword differs between MATRIX and EBOUNDS extensions.";
    SPreportError(InconsistentKeywordValues, msg);
  }
  if ( m_Instrument != eBounds.getInstrument() ) {
    string msg = "INSTRUME keyword differs between MATRIX and EBOUNDS extensions.";
    SPreportError(InconsistentKeywordValues, msg);
  }
  if ( m_Detector != eBounds.getDetector() ) {
    string msg = "DETNAM keyword differs between MATRIX and EBOUNDS extensions.";
    SPreportError(InconsistentKeywordValues, msg);
  }
  if ( m_Filter != eBounds.getFilter() ) {
    string msg = "FILTER keyword differs between MATRIX and EBOUNDS extensions.";
    SPreportError(InconsistentKeywordValues, msg);
  }

  // copy the relevant pieces from eBounds to this.

  m_EBDExtensionName = eBounds.getEBDExtensionName();
  Status = this->setChannelLowEnergy(eBounds.getChannelLowEnergy());
  Status = this->setChannelHighEnergy(eBounds.getChannelHighEnergy());

  return(Status);
}

// reading channel bounds extension from PHA file. this option allows multiple extensions in the same file

Integer rmf::readChannelBounds(const string filename, const Integer EBDnumber)
{
  bool verbosity = FITS::verboseMode();

  // Attempt to open filename as a FITS object

  unique_ptr<FITS> pInfile((FITS*)0);
  try {
    pInfile.reset(new FITS(filename,Read));
  } catch(...) {
    string msg = "Failed to open "+filename;
    SPreportError(NoSuchFile, msg);
    return(NoSuchFile);
  }

  // File exists so now look for the matrix extension. This can be identified
  // by the hduName being "EBOUNDS". If this doesn't work also check for the 
  // presence of the HDUCLAS1 and HDUCLAS2 keywords being set to "RESPONSE" 
  // and "EBOUNDS", respectively.

  string hduName("EBOUNDS");
  const vector<string> hduKeys;
  try {
    pInfile->read(hduName, Read, hduKeys, EBDnumber);
  } catch(...) {
    try {
      vector<string> searchKeys(2);
      vector<string> searchValues(2);
      searchKeys[0] = "HDUCLAS1";
      searchKeys[1] = "HDUCLAS2";
      searchValues[0] = "RESPONSE";
      searchValues[1] = "EBOUNDS";
      pInfile->read(searchKeys,searchValues, Read, hduKeys, EBDnumber);
    } catch(...) {
      string msg = "Failed to find energy bounds extension in "+filename;
      SPreportError(NoSuchFile, msg);
      return(NoSuchFile);
    }
  }

  // if we have got to this point the current extension must be the
  // one we need

  ExtHDU& ebd = pInfile->currentExtension();

  // read the version number in case we need the information

  string DefString = "UNKNOWN";

  string EBDVersion = SPreadKey(ebd, "HDUVERS", DefString);
  if ( EBDVersion == "UNKNOWN" ) {
    EBDVersion = SPreadKey(ebd, "HDUVERS2", DefString);
    if ( EBDVersion == "UNKNOWN" ) {
      EBDVersion = SPreadKey(ebd, "EBDVERSN", DefString);
      if ( EBDVersion == "UNKNOWN" ) EBDVersion = "1.1.0";
    }
  }

  // read the standard keywords and store in the object

  m_ChannelType = SPreadKey(ebd, "CHANTYPE", DefString);

  m_EBDExtensionName = SPreadKey(ebd, "EXTNAME", DefString);
  
  m_Telescope = SPreadKey(ebd, "TELESCOP", DefString);
  
  m_Instrument = SPreadKey(ebd, "INSTRUME", DefString);

  m_Detector = SPreadKey(ebd, "DETNAM", DefString);

  m_Filter = SPreadKey(ebd, "FILTER", DefString);

  // Read the E_MIN and E_MAX columns

  SPreadCol(ebd,"E_MIN",m_ChannelLowEnergy);
  SPreadCol(ebd,"E_MAX",m_ChannelHighEnergy);
  if ( m_ChannelLowEnergy.size() == 0 ) {
    string msg = "Failed to read any entries from E_MIN";
    SPreportError(NoEmin, msg);
    return(NoEmin);
  }
  if ( m_ChannelHighEnergy.size() == 0 ) {
    string msg = "Failed to read any entries from E_MAX";
    SPreportError(NoEmax, msg);
    return(NoEmax);
  }

  SPreadColUnits(ebd, "E_MIN", m_EnergyUnits);

  FITS::clearErrors();
  FITS::setVerboseMode(verbosity);

  return(OK);
}

// reading Matrix extension from PHA file.

Integer rmf::readMatrix(const string filename, const Integer RMFnumber)
{
  bool verbosity = FITS::verboseMode();

  // Attempt to open filename as a FITS object

  unique_ptr<FITS> pInfile((FITS*)0);
  try {
    pInfile.reset(new FITS(filename,Read));
  } catch(...) {
    string msg = "Failed to open "+filename;
    SPreportError(NoSuchFile, msg);
    return(NoSuchFile);
  }

  // File exists so now look for the matrix extension. This can be identified
  // by the hduName being "MATRIX" or "SPECRESP MATRIX". If neither of those
  // works also check for the presence of the HDUCLAS1 and HDUCLAS2 keywords 
  // being set to "RESPONSE" and "RSP_MATRIX", respectively.

  string hduName("MATRIX");
  const vector<string> hduKeys;
  try {
    pInfile->read(hduName, Read, hduKeys, RMFnumber);
  } catch(...) {
    try {
      hduName = "SPECRESP MATRIX";
      pInfile->read(hduName, Read, hduKeys, RMFnumber);
    } catch(...) {
      try {
	vector<string> searchKeys(2);
	vector<string> searchValues(2);
	searchKeys[0] = "HDUCLAS1";
	searchKeys[1] = "HDUCLAS2";
	searchValues[0] = "RESPONSE";
	searchValues[1] = "RSP_MATRIX";
	pInfile->read(searchKeys,searchValues, Read, hduKeys, RMFnumber);
      } catch(...) {
	string msg = "Failed to find response matrix extension in "+filename;
	SPreportError(NoSuchFile, msg);
	return(NoSuchFile);
      }
    }
  }

  // if we have got to this point the current extension must be the
  // one we need

  ExtHDU& rmf = pInfile->currentExtension();

  // read the version number in case we need the information

  string DefString = "UNKNOWN";

  string RMFVersion = SPreadKey(rmf, "HDUVERS", DefString);
  if ( RMFVersion == "UNKNOWN" ) {
    RMFVersion = SPreadKey(rmf, "HDUVERS2", DefString);
    if ( RMFVersion == "UNKNOWN" ) {
      RMFVersion = SPreadKey(rmf, "RMFVERSN", DefString);
      if ( RMFVersion == "UNKNOWN" ) RMFVersion = "1.3.0";
    }
  }

  // read the standard keywords and store in the object

  m_ChannelType = SPreadKey(rmf, "CHANTYPE", DefString);

  m_RMFExtensionName = SPreadKey(rmf, "EXTNAME", DefString);
  
  m_Telescope = SPreadKey(rmf, "TELESCOP", DefString);
  
  m_Instrument = SPreadKey(rmf, "INSTRUME", DefString);

  m_Detector = SPreadKey(rmf, "DETNAM", DefString);

  m_Filter = SPreadKey(rmf, "FILTER", DefString);

  m_RMFType = SPreadKey(rmf, "HDUCLAS3", DefString);

  m_AreaScaling = SPreadKey(rmf, "EFFAREA", (Real)1.0);

  m_ResponseThreshold = SPreadKey(rmf, "LO_THRES", (Real)0.0);

  Integer Nrows;
  Nrows = SPreadKey(rmf, "NAXIS2", (Integer)0);

  Integer NtotGroups;
  NtotGroups = SPreadKey(rmf, "NUMGRP", (Integer)0);
  
  Integer NtotElts;
  NtotElts = SPreadKey(rmf, "NUMELT", (Integer)0);
 
  // Read the ENERG_LO and ENERG_HI columns

  SPreadCol(rmf,"ENERG_LO",m_LowEnergy);
  SPreadCol(rmf,"ENERG_HI",m_HighEnergy);
  if ( m_LowEnergy.size() == 0 ) {
    string msg = "Failed to read any entries from the ENERG_LO column";
    SPreportError(NoEnergLo, msg);
    return(NoEnergLo);
  }
  if ( m_HighEnergy.size() == 0 ) {
    string msg = "Failed to read any entries from the ENERG_HI column";
    SPreportError(NoEnergHi, msg);
    return(NoEnergHi);
  }

  SPreadColUnits(rmf, "ENERG_LO", m_EnergyUnits);

  // Get the number of groups for each energy bin

  SPreadCol(rmf,"N_GRP",m_NumberGroups);
  if ( m_NumberGroups.size() == 0 ) {
    string msg = "Failed to read any entries from the N_GRP column";
    SPreportError(NoNgrp, msg);
    return(NoNgrp);
  }

  // Test for consistency between any value read from the NUMGRP keyword and the
  // sum of the N_GRP column.

  Integer Ntest(0);
  for (size_t i=0; i<(size_t)Nrows; i++) Ntest += m_NumberGroups[i];

  if ( NtotGroups != 0  && Ntest != NtotGroups ) {
    stringstream msg;
    msg << "The value of NUMGRP (" << NtotGroups << ") does not equal the sum of the N_GRP column (" << Ntest << "), using the latter.";
    SPreportError(InconsistentNumgrp, msg.str());
  }
  NtotGroups = Ntest;

  // Read the first channel and number of channels for each group

  SPreadVectorCol(rmf, "F_CHAN", m_FirstChannelGroup);
  SPreadVectorCol(rmf, "N_CHAN", m_NumberChannelsGroup);
  if ( m_FirstChannelGroup.size() == 0 ) {
    string msg = "Failed to read any entries from the F_CHAN column";
    SPreportError(NoFchan, msg);
    return(NoFchan);
  }
  if ( m_NumberChannelsGroup.size() == 0 ) {
    string msg = "Failed to read any entries from the N_CHAN column";
    SPreportError(NoNchan, msg);
    return(NoNchan);
  }

  // Sometimes an input response file will include channel group information
  // with first channel and number of channels 0 even if there are no channel 
  // groups for this energy. Tidy up this case by resizing the appropriate arrays
  // to zero.

  for (size_t i=0; i<m_NumberGroups.size(); i++) {
    if ( m_NumberGroups[i] == 0 ) {
      m_FirstChannelGroup[i].resize(0);
      m_NumberChannelsGroup[i].resize(0);
    }
  }

  // Check for a TLMIN for the F_CHAN column

  try {
    Integer ChannelIndex = rmf.column("F_CHAN").index();
    ostringstream KeyStream;
    KeyStream << "TLMIN" << ChannelIndex;
    rmf.readKey(KeyStream.str(),m_FirstChannel);
    m_FirstChannel = SPreadKey(rmf, KeyStream.str(), (Integer)1);
  } catch(Table::NoSuchColumn&) {
    m_FirstChannel = 1;
  } catch(HDU::NoSuchKeyword&) {
    m_FirstChannel = 1;
  }

  // Test for consistency between any value read from the NUMELT keyword and the
  // sum of the N_CHAN column.

  Ntest = 0;
  for (size_t i=0; i<(size_t)m_NumberChannelsGroup.size(); i++) {
    for (size_t j=0; j<(size_t)m_NumberChannelsGroup[i].size(); j++) {
      Ntest += m_NumberChannelsGroup[i][j];
    }
  }

  if ( NtotElts != 0 && Ntest != NtotElts ) {
    stringstream msg;
    msg << "The value of NUMELT (" << NtotElts << ") does not equal the sum of the N_CHAN column (" << Ntest << "), using the latter.";
    SPreportError(InconsistentNumelt, msg.str());
  }
  NtotElts = Ntest;

  // Read the matrix column

  SPreadVectorCol(rmf, "MATRIX", m_Matrix);
  if ( m_Matrix.size() == 0 ) {
    string msg = "Failed to read any entries from the MATRIX column";
    SPreportError(NoMatrix, msg);
    return(NoMatrix);
  }

  // We have to be a bit careful here in case the file was created using fixed
  // length vectors. In this case the total size of the elements matrix read
  // will be larger than NtotElts with zero padding to the right in each row.
  // Modify m_Matrix to remove the padding

  for (size_t i=0; i<(size_t)m_Matrix.size(); i++) {
    size_t NtoRead(0);
    for (size_t j=0; j<(size_t)m_NumberGroups[i]; j++) {
      NtoRead += m_NumberChannelsGroup[i][j];
    }
    if ( m_Matrix[i].size() != NtoRead ) {
      vector<Real> saveArray(m_Matrix[i]);
      m_Matrix[i].resize(NtoRead);
      for (size_t j=0; j<NtoRead; j++) {
	m_Matrix[i][j] = saveArray[j];
      }
    }
  }

  SPreadColUnits(rmf, "MATRIX", m_RMFUnits);

  // Read the optional order information

  try {
    FITS::setVerboseMode(false);
    SPreadVectorCol(rmf, "ORDER", m_OrderGroup);
  } catch (...) {
  }

  FITS::clearErrors();
  FITS::setVerboseMode(verbosity);

  return(OK);
}

// Load object from input information

Integer rmf::load(const vector<Integer>& inNumberGroups,
		  const vector<vector<Integer> >& inFirstChannelGroup,
		  const vector<vector<Integer> >& inNumberChannelsGroup,
		  const vector<vector<Integer> >& inOrderGroup,
		  const vector<Real>& inLowEnergy, const vector<Real>& inHighEnergy,
		  const vector<vector<Real> >& inMatrix, 
		  const vector<Real>& inChannelLowEnergy,
		  const vector<Real>& inChannelHighEnergy,
		  const Integer inFirstChannel, const Real inAreaScaling,
		  const Real inResponseThreshold, const map<string,string>& inKeys)
{
  setNumberGroups(inNumberGroups);
  setFirstChannelGroup(inFirstChannelGroup);
  setNumberChannelsGroup(inNumberChannelsGroup);
  setOrderGroup(inOrderGroup);
  setLowEnergy(inLowEnergy);
  setHighEnergy(inHighEnergy);
  setMatrix(inMatrix);
  setChannelLowEnergy(inChannelLowEnergy);
  setChannelHighEnergy(inChannelHighEnergy);

  setFirstChannel(inFirstChannel);
  setAreaScaling(inAreaScaling);
  setResponseThreshold(inResponseThreshold);

  try { setEnergyUnits(inKeys.at("EnergyUnits")); } catch(...) {};
  try { setRMFUnits(inKeys.at("RMFUnits")); } catch(...) {};
  try { setChannelType(inKeys.at("ChannelType")); } catch(...) {};
  try { setTelescope(inKeys.at("Telescope")); } catch(...) {};
  try { setInstrument(inKeys.at("Instrument")); } catch(...) {};
  try { setDetector(inKeys.at("Detector")); } catch(...) {};
  try { setFilter(inKeys.at("Filter")); } catch(...) {};
  try { setRMFType(inKeys.at("RMFType")); } catch(...) {};
  try { setRMFExtensionName(inKeys.at("RMFExtensionName")); } catch(...) {};
  try { setEBDExtensionName(inKeys.at("EBDExtensionName")); } catch(...) {};

  return(OK);
}

// initialize from an arf object. Copies members in common between arfs and rmfs

void rmf::initialize(const arf& a)
{
  m_EnergyUnits = a.getEnergyUnits();
  m_Telescope   = a.getTelescope();
  m_Instrument  = a.getInstrument();
  m_Detector    = a.getDetector();
  m_Filter      = a.getFilter();

  vector<Real> arfLowEnergy = a.getLowEnergy();
  vector<Real> arfHighEnergy = a.getHighEnergy();
  
  m_LowEnergy.resize(arfLowEnergy.size());
  for (size_t i=0; i<m_LowEnergy.size(); i++) m_LowEnergy[i] = arfLowEnergy[i];
  m_HighEnergy.resize(arfHighEnergy.size());
  for (size_t i=0; i<m_HighEnergy.size(); i++) m_HighEnergy[i] = arfHighEnergy[i];

  return;
}

void rmf::loadDiagonalResponse(const vector<Real>& eLow, const vector<Real>& eHigh,
			       const vector<Real>& rspVals, const Integer firstChan)
{
  m_FirstChannel = firstChan;

  size_t nBins = eLow.size();

  m_LowEnergy.resize(nBins);
  for (size_t i=0; i<nBins; i++) m_LowEnergy[i] = eLow[i];
  m_HighEnergy.resize(nBins);
  for (size_t i=0; i<nBins; i++) m_HighEnergy[i] = eLow[i];

  m_ChannelLowEnergy.resize(nBins);
  for (size_t i=0; i<nBins; i++) m_LowEnergy[i] = eLow[i];
  m_ChannelHighEnergy.resize(eLow.size());
  for (size_t i=0; i<nBins; i++) m_HighEnergy[i] = eLow[i];

  m_NumberGroups.resize(nBins);
  m_FirstChannelGroup.resize(nBins);
  m_NumberChannelsGroup.resize(nBins);
  m_Matrix.resize(nBins);

  for (size_t i=0; i<nBins; i++) {
    m_NumberGroups[i] = 1;
    m_FirstChannelGroup[i].resize(1);
    m_FirstChannelGroup[i][0] = i+m_FirstChannel;
    m_NumberChannelsGroup[i].resize(1);
    m_NumberChannelsGroup[i][0] = 1;
    m_Matrix[i].resize(1);
    m_Matrix[i][0] = rspVals[i];
  }

  return;
}


// Return the value for a particular channel, energy and order. Use special case of
// GratingOrder = -999 for ignore order.

Real rmf::ElementValue(const Integer ChannelNumber, const Integer EnergyBin,
		       const Integer GratingOrder) const
{                                          

  if ( ChannelNumber < m_FirstChannel || ChannelNumber >= m_FirstChannel+NumberChannels() || 
       EnergyBin < 0 || EnergyBin >= NumberEnergyBins() ) return 0.0;

  // loop round groups for this energy bin

  for(size_t i=0;i<(size_t)(m_NumberGroups[EnergyBin]);i++) {

    if ( ( m_OrderGroup.size() > 0 && m_OrderGroup[EnergyBin][i] == GratingOrder ) 
	 || GratingOrder == -999 ) {

      if( ChannelNumber >= m_FirstChannelGroup[EnergyBin][i] && 
	  ChannelNumber < m_FirstChannelGroup[EnergyBin][i]+m_NumberChannelsGroup[EnergyBin][i]) {

	return(m_Matrix[EnergyBin][ChannelNumber-m_FirstChannelGroup[EnergyBin][i]]);

      }

    }

  }

  return(0.0);

}

// Return vector of matrix values for a particular energy and grating order
// Use GratingOrder = -999 as a special case for no check.

vector<Real> rmf::RowValues(const Integer EnergyBin, const Integer GratingOrder) const
{

  vector<Real> values(NumberChannels(),0.0);

  // Loop round response groups for this energy

  size_t ielt(0);

  for (size_t igroup=0; igroup<(size_t)m_NumberGroups[EnergyBin]; igroup++) {

    // loop round elements in this group - adding them to the output array

    if ( ( m_OrderGroup.size() > 0 && m_OrderGroup[EnergyBin][igroup] == GratingOrder ) 
	 || GratingOrder == -999 ) {

      size_t ivec(m_FirstChannelGroup[EnergyBin][igroup]-m_FirstChannel);

      for (size_t j=0; j<(size_t)m_NumberChannelsGroup[EnergyBin][igroup]; j++) {
	values[ivec+j] += m_Matrix[EnergyBin][ielt+j];
      }

      ielt += m_NumberChannelsGroup[EnergyBin][igroup];

    }

  }

  return values;
}

// Return vector of randomly generated channel numbers for a particular energy and
// grating order. Use GratingOrder = -999 as special case to ignore grating information

vector<Integer> rmf::RandomChannels(const Real energy, const Integer NumberPhotons,
				    const vector<Real>& RandomNumber,
				    const Integer GratingOrder) const
{
  return this->RandomChannels(vector<Real>(1,energy), vector<Integer>(1,NumberPhotons),
			      vector<vector<Real> >(1,RandomNumber), GratingOrder);
}

// Return vector of randomly generated channel numbers for a set of energies and a
// grating order. Use GratingOrder = -999 as special case to ignore grating information

vector<Integer> rmf::RandomChannels(const vector<Real>& energy,
				    const vector<Integer>& NumberPhotons,
				    const vector<vector<Real> >& RandomNumber,
				    const Integer GratingOrder) const
{
  Integer NumberOut(0);
  for (size_t i=0; i<NumberPhotons.size(); i++) NumberOut += NumberPhotons[i];
  vector<Integer> channel(NumberOut);

  // initialize the output array to -1s in the event that either the input energy is
  // outside the response range or that the response does not sum to unity and events
  // can fall off the end of the channels.

  for (size_t i=0; i<(size_t)NumberOut; i++) channel[i] = -1;

  // loop round the energies

  size_t iout(0); 
  Real Emin = m_LowEnergy[0];
  Real Emax = m_HighEnergy[m_HighEnergy.size()-1];

  for (size_t i=0; i<energy.size(); i++) {

    // trap the case of the energy being outside the response range

    if ( energy[i] >= Emin && energy[i] <= Emax ) {

      // find the energy bin associated with the input energy 
      // - assumes the energies are in increasing order

      size_t lower = 0;
      size_t upper = m_HighEnergy.size()-1;
      size_t middle, energybin;
      while ( upper - lower > 1 ) {
	middle = (upper + lower)/2;
	if ( energy[i] < m_HighEnergy[middle] ) {
	  upper = middle;
	} else {
	  lower = middle;
	}
      }
      if ( energy[i] > m_HighEnergy[lower] ) {
	energybin = upper;
      } else {
	energybin = lower;
      }

      // generate an array of size channel each element of which is the integrated 
      // response up to and including that channel

      vector<Real> sumresponse(m_ChannelHighEnergy.size());
      for (size_t j=0; j<m_ChannelHighEnergy.size(); j++) sumresponse[j] = 0.0;

      sumresponse = this->RowValues(energybin, GratingOrder);

      for (size_t j=1; j<sumresponse.size(); j++) sumresponse[j] += sumresponse[j-1];

      // loop round the photons

      for (size_t j=0; j<(size_t)RandomNumber[i].size(); j++) {

	Real random = RandomNumber[i][j];

	// find the array element containing this random number. note that we do
	// not assume that the total response sums to 1 - if the random number
	// exceeds the total response then we assume that the event fell off the
	// end of the channel array and return a -1

	lower = 0;
	upper = m_ChannelHighEnergy.size() - 1;
	if ( random <= sumresponse[upper] ) {
	  while ( upper - lower > 1 ) {
	    middle = (upper + lower)/2;
	    if ( random < sumresponse[middle] ) {
	      upper = middle;
	    } else {
	      lower = middle;
	    }
	  }
	  if ( random > sumresponse[lower] ) {
	    channel[iout] = upper;
	  } else {
	    channel[iout] = lower;
	  }

	  // correct the channel number for the first channel number in use in 
	  // the response matrix

	  channel[iout] += m_FirstChannel;

	}
	iout++;

	// end loop over photons
      }

    } else {

      // if the current energy is outside the response range then increment iout by
      // the number of photons we would be simulating

      iout += NumberPhotons[i];

    }

    // end loop over energies

  }

  return channel;
}

// Display information about the RMF - return as a string

string rmf::disp() const
{
  ostringstream outstr;

  outstr <<"Response information : " << endl;

  outstr << "   FirstChannel        = " << m_FirstChannel << endl;
  outstr << "   AreaScaling         = " << m_AreaScaling << endl;
  outstr << "   ResponseThreshold   = " << m_ResponseThreshold << endl;
  outstr << "   ChannelType         = " << m_ChannelType << endl;
  outstr << "   Telescope           = " << m_Telescope << endl;
  outstr << "   Instrument          = " << m_Instrument << endl;
  outstr << "   Detector            = " << m_Detector << endl;
  outstr << "   Filter              = " << m_Filter << endl;
  outstr << "   RMFType             = " << m_RMFType << endl;
  outstr << "   RMFExtensionName    = " << m_RMFExtensionName << endl;
  outstr << "   EBDExtensionName    = " << m_EBDExtensionName << endl;

  outstr << "   EnergyUnits         = " << m_EnergyUnits << endl;
  outstr << "   RMFUnits            = " << m_RMFUnits << endl;

  outstr << "   NumberChannels      = " << NumberChannels() << endl;
  outstr << "   NumberEnergyBins    = " << NumberEnergyBins() << endl;
  outstr << "   NumberTotalGroups   = " << NumberTotalGroups() << endl;
  outstr << "   NumberTotalElements = " << NumberTotalElements() << endl;

  return outstr.str();
}

// Clear information from the response

void rmf::clear()
{
  m_FirstChannel = 0;

  this->clearMatrix();

  m_ChannelLowEnergy.clear();
  m_ChannelHighEnergy.clear();

  m_AreaScaling = 0.0;
  m_ResponseThreshold = 0.0;

  m_EnergyUnits = " ";
  m_RMFUnits = " ";

  m_ChannelType = " ";
  m_Telescope = " ";
  m_Instrument = " ";
  m_Detector = " ";
  m_Filter = " ";
  m_RMFType = " ";
  m_RMFExtensionName = " ";
  m_EBDExtensionName = " ";

  return;
}

// Clear only the matrix from the response

void rmf::clearMatrix()
{
  m_NumberGroups.clear();
  m_FirstChannelGroup.clear();
  m_NumberChannelsGroup.clear();
  m_OrderGroup.clear();
  m_LowEnergy.clear();
  m_HighEnergy.clear();
  m_Matrix.clear();

  return;
}

// Check completeness and consistency of information in the rmf
  // if there is a problem then return diagnostic in string

string rmf::check() const
{
  ostringstream outstr;

  // check for presence of any data

  if ( m_Matrix.size() == 0 ) {
    outstr << "Matrix has no data" << endl;
  }

  // check size consistency between arrays - channels

  if ( m_ChannelLowEnergy.size() != m_ChannelHighEnergy.size() ) {
    outstr << "ChannelLowEnergy size (" << m_ChannelLowEnergy.size() 
	 << ") differs from ChannelHighEnergy size (" << m_ChannelHighEnergy.size() 
	 << ")" << endl;
  }

  // energies

  if ( m_LowEnergy.size() != m_HighEnergy.size() ) {
    outstr << "LowEnergy size (" << m_LowEnergy.size() 
	 << ") differs from HighEnergy size (" << m_HighEnergy.size() 
	 << ")" << endl;
  }

  if ( m_LowEnergy.size() != m_NumberGroups.size() ) {
    outstr << "LowEnergy size (" << m_LowEnergy.size() 
	 << ") differs from NumberGroups size (" << m_NumberGroups.size() 
	 << ")" << endl;
  }

  // groups

  size_t tot1(0), tot2(0), tot3(0);
  for (size_t i=0; i<m_FirstChannelGroup.size(); i++) tot1 += m_FirstChannelGroup[i].size();
  for (size_t i=0; i<m_NumberChannelsGroup.size(); i++) tot2 += m_NumberChannelsGroup[i].size();
  for (size_t i=0; i<m_OrderGroup.size(); i++) tot3 += m_OrderGroup[i].size();

  if ( tot1 != tot2 ) {
    outstr << "The total number of FirstChannelGroup entries (" << tot1 
	 << ") differs from that for NumberChannelsGroup (" << tot2 << ")" << endl;
  }

  if ( tot1 != tot3 && m_OrderGroup.size() > 0 ) {
    outstr << "The total number of FirstChannelGroup entries (" << tot1
	 << ") differs from that for OrderGroup (" << tot3 << ")" << endl;
  }

  // check that arrays have sensible values

  if ( tot2  == 0 ) {
    outstr << "No groups have any channels - something went very wrong" << endl;
  } else {
    for (size_t i=0; i<m_NumberGroups.size(); i++) {
      if ( m_NumberGroups[i] < 0 ) {
	outstr << "NumberGroups has invalid value (" << m_NumberGroups[i] 
	       << ") for energy bin " << i << endl;
      }
    }
  }

  for (size_t i=0; i<m_FirstChannelGroup.size(); i++) {
    for (size_t j=0; j<m_FirstChannelGroup[i].size(); j++) {
      if ( m_FirstChannelGroup[i][j] < m_FirstChannel || 
	   m_FirstChannelGroup[i][j] >= (Integer)m_ChannelLowEnergy.size() + m_FirstChannel ) {
	outstr << "FirstChannelGroup has invalid value (" << m_FirstChannelGroup[i][j] 
	       << ") for the " << j << " group of energy bin " << i 
	       <<  ". Should be >= " << m_FirstChannel << " and < " 
	       << m_ChannelLowEnergy.size() <<endl;
      }
      if ( m_NumberChannelsGroup[i][j] < 0 || m_NumberChannelsGroup[i][j] > (Integer)m_ChannelLowEnergy.size() ) {
	outstr << "NumberChannelsGroup has invalid value (" << m_NumberChannelsGroup[i][j] 
	       << ") for the " << j <<  " group of energy bin " << i
	       << ". Should be >= 0 and <= " 
	       << m_ChannelLowEnergy.size() <<endl;
      }
    }      
  }
  return outstr.str();
}

// Normalize the rmf so it sums to 1.0 for each energy bin
// Returns the normalization factors.

vector<Real> rmf::normalize()
{

  vector<Real> factors(NumberEnergyBins(),0.0);
  
  // Loop over energies

  for (size_t ie=0; ie<(size_t)NumberEnergyBins(); ie++) {

    // sum up the response in this energy

    for (size_t i=0; i<m_Matrix[ie].size(); i++) {
      factors[ie] += m_Matrix[ie][i];
    }

    // divide through by the summed response

    for (size_t i=0; i<m_Matrix[ie].size(); i++) {
      m_Matrix[ie][i] /= factors[ie];
    }

  }

  return factors;
}

void rmf::changeFirstChannel(const Integer first)
{
  // Find the number that will be added to go from the
  // current first channel number to that requested
  Integer change = first - m_FirstChannel;

  // reset the member specifying the first channel
  m_FirstChannel = first;

  // loop through the FirstChannelGroup array resetting
  for (size_t ie=0; ie<m_FirstChannelGroup.size(); ie++) {
    for (size_t igrp=0; igrp<m_FirstChannelGroup[ie].size(); igrp++) {
      m_FirstChannelGroup[ie][igrp] += change;
    }
  }

  return;
}

void rmf::compress(const Real threshold)
{

  // Reset response threshold

  this->setResponseThreshold(threshold);

  // Temporary array for the response for a given energy

  vector<Real> Response(m_ChannelLowEnergy.size());

  // Loop over energies

  for (size_t i=0; i<(size_t)m_LowEnergy.size(); i++) {

    // expand response matrix into a channel array for this energy

    Response = this->RowValues(i);

    // and replace in the response using the new threshold

    this->substituteRow(i, Response);

    // end loop over energies

  }

  return;
}

void rmf::uncompress()
{
  // can do this just by calling compress with a threshold of zero.

  this->compress(0.0);
  return;
}

// Compress in channel space

Integer rmf::rebinChannels(grouping& GroupInfo)
{

  // check for consistency between grouping and number of channels

  if ( GroupInfo.size() != NumberChannels() ) {
    stringstream msg;
    msg << "Number of channels (" << NumberChannels() << ") differs from size of grouping array (" << GroupInfo.size() << ").";
    SPreportError(InconsistentGrouping, msg.str());
    return(InconsistentGrouping);
  }

  // Temporary array for the response for a given energy

  vector<Real> Response(NumberChannels());

  // Loop over energies

  for (size_t i=0; i<(size_t)m_LowEnergy.size(); i++) {

    // expand response matrix into a channel array for this energy

    Response = this->RowValues(i);

    // bin up the response for the energy

    vector<Real> binResponse;
    GroupBin(Response, SumMode, GroupInfo, binResponse);

    // update the compressed response arrays for this energy bin

    this->substituteRow(i, binResponse);

    // end loop over energies

  }

  // now rebin the channel boundary arrays
  // the 2 arrays may be in ascending or descending order
  
  const size_t nEboundsEngs = m_ChannelLowEnergy.size();
  const bool isDescending = (nEboundsEngs && 
      (m_ChannelLowEnergy[nEboundsEngs-1] < m_ChannelLowEnergy[0]));
  
  vector<Real> temp;
  if (isDescending)
     GroupBin(m_ChannelLowEnergy, LastEltMode, GroupInfo, temp);
  else
     GroupBin(m_ChannelLowEnergy, FirstEltMode, GroupInfo, temp);
  m_ChannelLowEnergy.resize(temp.size());
  for (size_t i=0; i<m_ChannelLowEnergy.size(); i++) {
    m_ChannelLowEnergy[i] = temp[i];
  }
  if (isDescending)
     GroupBin(m_ChannelHighEnergy, FirstEltMode, GroupInfo, temp);
  else
     GroupBin(m_ChannelHighEnergy, LastEltMode, GroupInfo, temp);
  m_ChannelHighEnergy.resize(temp.size());
  for (size_t i=0; i<m_ChannelHighEnergy.size(); i++) {
    m_ChannelHighEnergy[i] = temp[i];
  }

  return(OK);
}

// Compress in energy space

Integer rmf::rebinEnergies(grouping& GroupInfo)
{

  // check for consistency between grouping and number of energy bins

  if ( GroupInfo.size() != NumberEnergyBins() ) {
    stringstream msg;
    msg << "Number of energy bins (" << NumberEnergyBins() << ") differs from size of grouping array (" << GroupInfo.size() << ").";
    SPreportError(InconsistentGrouping, msg.str());
    return(InconsistentGrouping);
  }

  // Set up temporary object to store the output rmf and set its threshold and first channel

  rmf work;
  work.setResponseThreshold(m_ResponseThreshold);
  work.setFirstChannel(m_FirstChannel);

  // Temporary array for the response for a given energy

  vector<Real> Response(NumberChannels(),0.0);
  Real eLow(m_LowEnergy[0]);
  Real eHigh(m_HighEnergy[0]);
  size_t nbins(0);

  // Loop over energies

  for (size_t i=0; i<(size_t)NumberEnergyBins(); i++) {

    // reset the Response array if this is the start of a new energy bin
    // and update the compressed response arrays if this is not the first bin

    if ( GroupInfo.newBin(i) ) {
      if ( i != 0 ) {
	if ( m_RMFType == "REDIST" && nbins > 0 ) {
	  for (size_t j=0; j<Response.size(); j++) Response[j] /= nbins;
	}
	work.addRow(Response, eLow, eHigh);
	nbins = 0;
      }
      for (size_t j=0; j<Response.size(); j++) Response[j] = 0.0;
      eLow = m_LowEnergy[i];
    }

    eHigh = m_HighEnergy[i];
    
    // expand response matrix into a channel array for this energy and accumulate

    vector<Real> RespArray(NumberChannels());
    RespArray = this->RowValues(i);
    for (size_t j=0; j<Response.size(); j++) Response[j] += RespArray[j];
    nbins++;

  }

  // update the compressed response arrays with the final row.

  if ( m_RMFType == "REDIST" && nbins > 0 ) {
    for (size_t j=0; j<Response.size(); j++) Response[j] /= nbins;
  }
  work.addRow(Response, eLow, eHigh);

  // copy new response arrays into current response

  setNumberGroups(work.getNumberGroups());
  setFirstChannelGroup(work.getFirstChannelGroup());
  setNumberChannelsGroup(work.getNumberChannelsGroup());
  setMatrix(work.getMatrix());
  setLowEnergy(work.getLowEnergy());
  setHighEnergy(work.getHighEnergy());

  return(OK);
}

// Shift channels up or down.

Integer rmf::shiftChannels(const Integer Start, const Integer End, const Real Shift, const Real Factor, bool useEnergyBounds)
{
  return this->shiftChannels(vector<Integer>(1,Start), vector<Integer>(1,End),
			     vector<Real>(1,Shift), vector<Real>(1,Factor),
			     useEnergyBounds);
}

Integer rmf::shiftChannels(const vector<Integer>& vStart, const vector<Integer>& vEnd, const vector<Real>& vShift, const vector<Real>& vFactor, bool useEnergyBounds)
{

  size_t Nchan(NumberChannels());
  size_t Nener(NumberEnergyBins());

  // First set up vectors describing how to make the new matrix
  // fromChannel[i] is the list of channels from which the new channel i is calculated
  // fromFraction[i] is the list of fractions corresponding to fromChannel.

  vector<vector<size_t> > fromChannel(Nchan);
  vector<vector<Real> > fromFraction(Nchan);

  if ( useEnergyBounds ) {
    SPcalcShift(m_ChannelLowEnergy, m_ChannelHighEnergy, vStart, vEnd, vShift, vFactor, 
		fromChannel, fromFraction);
  } else {
    // Shift is in terms of channel number so define the Low and High as the channel 
    //number -/+ 0.5
    vector<Real> Low(Nchan);
    vector<Real> High(Nchan);
    for (size_t i=0; i<Nchan; i++) {
      Low[i] = i + m_FirstChannel - 0.5;
      High[i] = i + m_FirstChannel + 0.5;
    }
    SPcalcShift(Low, High, vStart, vEnd, vShift, vFactor, fromChannel, fromFraction);
  }

  // Loop over energies accumulating the new response

  for (size_t iEnergyBin=0; iEnergyBin<Nener; iEnergyBin++) {

    vector<Real> Response(Nchan);
    vector<Real> OutResponse(Nchan, 0.0);
    
    // expand response matrix into a channel array for this energy

    Response = this->RowValues(iEnergyBin);

    // Construct the output response using the information in fromChannel
    // and fromFraction

    for ( size_t iChan=0; iChan<Nchan; iChan++) {

      for (size_t j=0; j<fromChannel[iChan].size(); j++) {
	OutResponse[iChan] += fromFraction[iChan][j]*Response[fromChannel[iChan][j]];
      }

    }

    // update the compressed response arrays for this energy bin

    this->substituteRow(iEnergyBin, OutResponse);

    // end loop over energies

  }

  return(OK);
}

// Shift energies up or down.

Integer rmf::shiftEnergies(const Integer Start, const Integer End, const Real Shift, const Real Factor)
{
  return this->shiftEnergies(vector<Integer>(1,Start), vector<Integer>(1,End),
			     vector<Real>(1,Shift), vector<Real>(1,Factor));
}

Integer rmf::shiftEnergies(const vector<Integer>& vStart, const vector<Integer>& vEnd, const vector<Real>& vShift, 
			   const vector<Real>& vFactor)
{

  // Note that Start and End are zero-based

  size_t Nchan(NumberChannels());
  size_t Nener(NumberEnergyBins());

  // Set up vectors describing how to make the new matrix
  // fromRow[i] is the list of rows contributing to row i
  // and Fraction[i] is the fractional contribution for each row

  vector<vector<size_t> > fromRow(Nener);
  vector<vector<Real> > fromFraction(Nener);

  SPcalcShift(m_LowEnergy, m_HighEnergy, vStart, vEnd, vShift, vFactor, fromRow, fromFraction);
 
  // Set up temporary object to store the output rmf and set its threshold 
  // and first channel

  rmf work;
  work.setResponseThreshold(m_ResponseThreshold);
  work.setFirstChannel(m_FirstChannel);

  // Loop over energies accumulating the new response

  for (size_t iEnergyBin=0; iEnergyBin<Nener; iEnergyBin++) {

    vector<Real> newRow(Nchan,0.0);

    // Construct the row from the contributions set in fromRow and fromFraction

    for (size_t i=0; i<fromRow[iEnergyBin].size(); i++) {

      vector<Real> oldRow(Nchan);
      size_t oldRowIndex = fromRow[iEnergyBin][i];
      oldRow = this->RowValues(oldRowIndex);
      Real frac = fromFraction[iEnergyBin][i];
      for (size_t k=0; k<Nchan; k++) newRow[k] += frac*oldRow[k];

    }

    work.addRow(newRow, m_LowEnergy[iEnergyBin], m_HighEnergy[iEnergyBin]);

  }

  // copy new response arrays into current response

  setNumberGroups(work.getNumberGroups());
  setFirstChannelGroup(work.getFirstChannelGroup());
  setNumberChannelsGroup(work.getNumberChannelsGroup());
  setMatrix(work.getMatrix());

  return(OK);
}

// Multiply by a vector which may not have the same energy binning as the response

Integer rmf::interpolateAndMultiply(const vector<Real>& inputEnergies, 
			       const vector<Real>& inputFactors)
{

  size_t nInputEnergies(inputEnergies.size());

  // loop over response energies

  for (size_t i=0; i<m_LowEnergy.size(); i++) {

    vector<Real> cx(2);
    cx[0] = m_LowEnergy[i];
    cx[1] = m_HighEnergy[i];

    // the multiplicative factor will be placed in factor

    Real factor;

    // If energy is above or below all input energy data then put factor equal to 0

    if ( cx[0] > inputEnergies[nInputEnergies-1] || cx[1] < inputEnergies[0] ) {

      factor = 0;

    } else {

      // set ibegin to the last input energy below the bottom of the bin

      size_t ibegin = 0;
      while ( inputEnergies[ibegin+1] < cx[0] ) ibegin++;

      // set iend to the first input energy above the top of the bin

      size_t iend = ibegin + 1;
      if ( iend >= nInputEnergies ) iend--;
      while ( inputEnergies[iend] < cx[1] && iend < nInputEnergies-1 ) iend++;

      // calculate number of tabulated values for this bin.

      size_t nValues = iend - ibegin + 1;

      // calculate interpolated values at top and bottom of bin. Catch cases
      // of the first or last input energies being in the bin. 

      vector<Real> cy(nValues);

      if ( cx[0] > inputEnergies[ibegin] ) {
	cy[0] = inputFactors[ibegin] + (inputFactors[ibegin+1]-inputFactors[ibegin])
	  *(cx[0]-inputEnergies[ibegin])/(inputEnergies[ibegin+1]-inputEnergies[ibegin]);
      } else {
	cy[0] = inputFactors[ibegin];
      }
      if ( cx[1] > inputEnergies[iend] ) {
	cy[nValues-1] = inputFactors[iend-1] + (inputFactors[iend]-inputFactors[iend-1])
	  *(cx[1]-inputEnergies[iend-1])/(inputEnergies[iend]-inputEnergies[iend-1]);
      } else {
	cy[nValues-1] = inputFactors[iend];
      }

      // if no input energies in current bin then factor is mean of these two

      if ( nValues <= 2 ) {

	factor = 0.5*(cy[0]+cy[1]);

      } else {

	// otherwise factor is energy-weighted mean

	cx.resize(nValues);
	cx[0] = m_LowEnergy[i];
	cx[nValues-1] = m_HighEnergy[i];

	for (size_t k=1; k<nValues-1; k++) {
	  cx[k] = inputEnergies[ibegin+k];
	  cy[k] = inputFactors[ibegin+k];
	}
	factor = 0.0;
	for (size_t k=1; k<nValues; k++) {
	  factor += 0.5*(cy[k]+cy[k-1])*(cx[k]-cx[k-1]);
	}
	factor /= (m_HighEnergy[i]-m_LowEnergy[i]);

      }

    }

    // multiply response for this energy by the factor

    for (size_t j=0; j<m_Matrix[i].size(); j++ ) m_Matrix[i][j] *= factor;

  }

  return(OK);

}
// Write response matrix and channel bounds extensions.

Integer rmf::write(const string filename, const string sourceFilename) const
{
  Integer Status(OK);

  // if we have a sourceFilename use that to set up the primary HDU for the
  // output rmf file otherwise just create a default primary HDU. Note that we
  // perform this entire operation within { } to clean up the objects before
  // calling the routines which make the EBOUNDS and MATRIX extensions.

  bool haveSourceFile(false);
  if ( sourceFilename.length() != 0 ) haveSourceFile = true;

  {
    std::unique_ptr<FITS> pSourceFits((FITS*)0);
    if ( haveSourceFile ) {
      try {
	pSourceFits.reset( new FITS(sourceFilename, Read, false) );
      } catch (...) {
	string msg = "Failed to read "+sourceFilename;
	SPreportError(NoSuchFile, msg);
	return(CannotCreate);
      }
    }

    // Create a new FITS file instance for the output file and try to create it

    std::unique_ptr<FITS> pFits((FITS*)0);
      
    try {
      if ( haveSourceFile ) {
	pFits.reset( new FITS(filename,*pSourceFits) );
      } else {
	pFits.reset( new FITS(filename, Write) );
      }
    } catch (FITS::CantCreate) {
      string msg = "Failed to create "+filename+" for response.";
      SPreportError(CannotCreate, msg);
      return(CannotCreate);       
    }

  }

  // these routines write the Matrix and ChannelBounds extensions.

  Status = this->writeChannelBounds(filename);
  if ( Status != OK ) return(Status);

  Status = this->writeMatrix(filename);
  if ( Status != OK ) return(Status);

  // Now if we had a source file input copy an extra keywords and extensions
  // into the output file. Slight subtlety here is that we need to check
  // whether the source file uses SPECRESP MATRIX or MATRIX.
  if ( haveSourceFile ) {
    Integer status = SPcopyKeys(sourceFilename, filename, "EBOUNDS");
    SPcopyCols(sourceFilename, filename, "EBOUNDS");
    status = SPcopyKeys(sourceFilename, filename, "MATRIX");
    if ( status ) {
      status = 0;
      status = SPcopyKeys(sourceFilename, filename, "SPECRESP MATRIX", "MATRIX");
      SPcopyCols(sourceFilename, filename, "SPECRESP MATRIX", "MATRIX");
    } else {
      SPcopyCols(sourceFilename, filename, "MATRIX");
    }
    status = SPcopyHDUs(sourceFilename, filename);
    return(status);
  }

  return(Status);
}

// Write response matrix extension. If file already exists appends.

Integer rmf::writeMatrix(const string filename) const
{
  string Blank = " ";

  vector<string> ttype;
  vector<string> tform;
  vector<string> tunit;

  // Create a new FITS file instance  

  std::unique_ptr<FITS> pFits((FITS*)0);
      
  try {                
    pFits.reset( new FITS(filename,Write) );
  } catch (FITS::CantCreate) {
    string msg = "Failed to create "+filename+" for MATRIX extension";
    SPreportError(CannotCreateMatrixExt, msg);
    return(CannotCreateMatrixExt);
  }

  // calculate the maximum number of groups and elements per row

  Integer Nrows = m_NumberGroups.size();
  
  Integer MaxGroups=0;
  Integer MaxElts=0;
  for (size_t i=0; i<m_NumberGroups.size(); i++) {
    if ( m_NumberGroups[i] > MaxGroups ) MaxGroups = m_NumberGroups[i];
    Integer NumElts(m_Matrix[i].size());
    if ( NumElts > MaxElts ) MaxElts = NumElts;
  }

  // set up the column descriptors for those attributes which need to be 
  // written as columns

  ttype.push_back("ENERG_LO");
  tform.push_back("D");
  tunit.push_back(m_EnergyUnits);

  ttype.push_back("ENERG_HI");
  tform.push_back("D");
  tunit.push_back(m_EnergyUnits);

  ttype.push_back("N_GRP");
  if ( MaxGroups > 32768 ) {
    tform.push_back("J");
  } else {
    tform.push_back("I");
  }
  tunit.push_back(" ");

  stringstream RepeatStream;
  RepeatStream << MaxGroups;
  string Repeat(RepeatStream.str());

  ttype.push_back("F_CHAN");
  if ( SPneedVecCol(m_FirstChannelGroup) ) {
    tform.push_back("PJ("+Repeat+")");
  } else {
    tform.push_back("J");
  }
  tunit.push_back(" ");

  ttype.push_back("N_CHAN");
  if ( SPneedVecCol(m_NumberChannelsGroup) ) {
    tform.push_back("PJ("+Repeat+")");
  } else {
    tform.push_back("J");
  }
  tunit.push_back(" ");

  bool isvector;
  if ( SPneedCol(m_OrderGroup, isvector) ) {
    ttype.push_back("ORDER");
    if ( isvector ) {
      tform.push_back("PJ("+Repeat+")");
    } else {
      tform.push_back("J");
    }
    tunit.push_back(" ");
  }

  RepeatStream.str("");
  RepeatStream << MaxElts;
  Repeat = RepeatStream.str();

  ttype.push_back("MATRIX");
  if ( SPneedVecCol(m_Matrix) ) {
    tform.push_back("PE("+Repeat+")");
  } else {
    tform.push_back("E");
  }
  tunit.push_back(m_RMFUnits);

  // Create the new extension. First check for existing extensions to see whether
  // we need to give a version number.

  ExtMapConstIt itLow = pFits->extension().lower_bound("MATRIX");
  ExtMapConstIt itHigh = pFits->extension().upper_bound("MATRIX");

  Integer version(0);
  while (itLow != itHigh) {
    if (itLow->second->version() > version) version = itLow->second->version();
    itLow++;
  }

  version++;
  Table* prmf = pFits->addTable("MATRIX",Nrows,ttype,tform,tunit,BinaryTbl,version);
  Table& rmf = *prmf;

  // Write the standard keywords
  
  SPwriteKey(rmf, "HDUCLASS", (string)"OGIP", Blank);
    
  SPwriteKey(rmf, "HDUCLAS1", (string)"RESPONSE", Blank);

  SPwriteKey(rmf, "HDUCLAS2", (string)"RSP_MATRIX", Blank);
    
  SPwriteKey(rmf, "HDUCLAS3", m_RMFType, Blank);
    
  SPwriteKey(rmf, "CHANTYPE", m_ChannelType, "Channel type");

  SPwriteKey(rmf, "HDUVERS", (string)"1.3.0", "OGIP version number");

  SPwriteKey(rmf, "TELESCOP", m_Telescope, Blank);

  SPwriteKey(rmf, "INSTRUME", m_Instrument, Blank);

  SPwriteKey(rmf, "DETNAM", m_Detector, Blank);

  SPwriteKey(rmf, "FILTER", m_Filter, Blank);

  SPwriteKey(rmf, "EFFAREA", m_AreaScaling, Blank);

  SPwriteKey(rmf, "LO_THRES", m_ResponseThreshold, Blank);

  SPwriteKey(rmf, "DETCHANS", NumberChannels(), "Number of channels in rmf");

  SPwriteKey(rmf, "NUMGRP", NumberTotalGroups(), "Total number of response groups");

  SPwriteKey(rmf, "NUMELT", NumberTotalElements(), "Total number of response elements");

  SPwriteKey(rmf, "TLMIN4", m_FirstChannel, "First channel number");

  // Write the arrays - if an array is of size 1 or all the same value 
  // it will be written as a keyword

  try {

    SPwriteCol(rmf, "ENERG_LO", m_LowEnergy, true);
    SPwriteCol(rmf, "ENERG_HI", m_HighEnergy, true);

    SPwriteCol(rmf, "N_GRP", m_NumberGroups, true);

    SPwriteVectorCol(rmf, "F_CHAN", m_FirstChannelGroup, true);
    SPwriteVectorCol(rmf, "N_CHAN", m_NumberChannelsGroup, true);

    SPwriteVectorCol(rmf, "MATRIX", m_Matrix, true);

    if ( m_OrderGroup.size() > 0 ) {
      SPwriteVectorCol(rmf, "ORDER", m_OrderGroup);
    }

  } catch(...) {

    string msg = "Failed to write MATRIX data to "+filename;
    SPreportError(CannotWriteMatrix, msg);
    return(CannotWriteMatrix);

  }

  // Write/update checksum keywords

  rmf.writeChecksum();

  return(OK);
}


// Write channel bounds extension. If file already exists appends.

Integer rmf::writeChannelBounds(const string filename) const
{

  string Blank = " ";

  vector<string> ttype;
  vector<string> tform;
  vector<string> tunit;

  // Create a new FITS file instance  

  std::unique_ptr<FITS> pFits((FITS*)0);
      
  try {                
    pFits.reset( new FITS(filename,Write) );
  } catch (FITS::CantCreate) {
    string msg = "Failed to create "+filename+" for EBOUNDS extension";
    SPreportError(CannotCreateEboundsExt, msg);
    return(CannotCreateEboundsExt);
  }

  // set up the column descriptors for those attributes which need to be 
  // written as columns

  ttype.push_back("CHANNEL");
  if ( m_ChannelLowEnergy.size() > 32768 ) {
    tform.push_back("J");
  } else {
    tform.push_back("I");
  }
  tunit.push_back(" ");

  ttype.push_back("E_MIN");
  tform.push_back("D");
  tunit.push_back(m_EnergyUnits);

  ttype.push_back("E_MAX");
  tform.push_back("D");
  tunit.push_back(m_EnergyUnits);

  // Create the new extension

  Table* pebd = pFits->addTable("EBOUNDS",m_ChannelLowEnergy.size(),ttype,tform,tunit);
  Table& ebd = *pebd;

  // Write the standard keywords
  
  SPwriteKey(ebd, "HDUCLASS", (string)"OGIP", Blank);
    
  SPwriteKey(ebd, "HDUCLAS1", (string)"RESPONSE", Blank);

  SPwriteKey(ebd, "HDUCLAS2", (string)"EBOUNDS", Blank);
    
  SPwriteKey(ebd, "CHANTYPE", m_ChannelType, "Channel type");

  SPwriteKey(ebd, "HDUVERS", (string)"1.1.0", "OGIP version number");

  SPwriteKey(ebd, "TELESCOP", m_Telescope, Blank);

  SPwriteKey(ebd, "INSTRUME", m_Instrument, Blank);

  SPwriteKey(ebd, "DETNAM", m_Detector, Blank);

  SPwriteKey(ebd, "FILTER", m_Filter, Blank);

  SPwriteKey(ebd, "DETCHANS", NumberChannels(), "Number of channels in ebd");

  // Generate and write the COLUMN array

  vector<Real> Channel(NumberChannels());
  for (size_t i=0; i<(size_t)NumberChannels(); i++) Channel[i] = i + m_FirstChannel;
  SPwriteCol(ebd, "CHANNEL", Channel, true);

  // Write the E_MIN and E_MAX arrays - if an array is of size 1 or all the 
  // same value it will be written as a keyword

  SPwriteCol(ebd, "E_MIN", m_ChannelLowEnergy, true);
  SPwriteCol(ebd, "E_MAX", m_ChannelHighEnergy, true);

  // Write/update checksum keywords

  ebd.writeChecksum();

  return(OK);
}


// Merge arf and rmf

rmf& rmf::operator*=(const arf& a)
{
  // check that the arf and rmf are compatible
  // if not just return the current rmf

  Integer status = checkCompatibility(a);
  if ( status != OK ) {
    SPreportError(status, "Failure when checking an RMF and an ARF for compatibility - they will not be multiplied.");
    return *this;
  }

  // loop round energy bins multiplying appropriate elements of the rmf by
  // the effective area for this energy from the ARF

  vector<Real> effarea = a.getEffArea();
  for ( size_t i=0; i < (size_t)m_LowEnergy.size(); i++ ) {
    for (size_t j=0; j < m_Matrix[i].size(); j++ ) m_Matrix[i][j] *= effarea[i];
  }

  return *this;
}

// Multiply by a constant factor

rmf& rmf::operator*=(const Real& f)
{
  for (size_t i=0; i<m_Matrix.size(); i++) {
    for (size_t j=0; j<m_Matrix[i].size(); j++) {
      m_Matrix[i][j] *= f;
    }
  }
  return *this;
}


rmf& rmf::operator+=(const rmf& r)
{

 // check that the two rmfs are compatible
 // if not just return the current rmf

  Integer status = checkCompatibility(r);
  if ( status != OK ) {
    SPreportError(status, "Failure when checking two RMFs for compatibility - they will not be summed.");
    return *this;
  }

  // temporary arrays for the response for each energy

  vector<Real> Response1(m_ChannelLowEnergy.size());
  vector<Real> Response2(m_ChannelLowEnergy.size());

  // loop round energy bins summing appropriate elements of the rmf

  for ( size_t i=0; i < (size_t)m_LowEnergy.size(); i++ ) {

    // expand both response matrices into a channel array for this energy

    Response1 = this->RowValues(i);
    Response2 = r.RowValues(i);

    // sum the two responses for this energy bin

    for (size_t j=0; j<Response1.size(); j++) {
      Response1[j] += Response2[j];
    }

    // update the compressed response arrays for this energy bin

    this->substituteRow(i, Response1);

    // end loop over energies

  }

  return *this;
}


Integer rmf::checkCompatibility(const rmf& r) const
{

  // check that the two rmf energy binnings are compatible

  vector<Real> rmfLowEnergy = r.getLowEnergy();
  vector<Real> rmfHighEnergy = r.getHighEnergy();

  if ( m_LowEnergy.size() != rmfLowEnergy.size() ) return InconsistentEnergies;
  for ( size_t i=0; i < (size_t)m_LowEnergy.size(); i++ ) {
    if ( fabs(m_LowEnergy[i]-rmfLowEnergy[i])/m_LowEnergy[i] > FUZZY ) return InconsistentEnergies;
    if ( fabs(m_HighEnergy[i]-rmfHighEnergy[i])/m_HighEnergy[i] > FUZZY ) return InconsistentEnergies;
  }

 // check that the two rmf channel binnings are compatible

  vector<Real> rmfChannelLowEnergy = r.getChannelLowEnergy();
  vector<Real> rmfChannelHighEnergy = r.getChannelHighEnergy();


  if ( m_ChannelLowEnergy.size() != rmfChannelLowEnergy.size() ) return InconsistentChannels;
  for ( size_t i=0; i < (size_t)m_ChannelLowEnergy.size(); i++ ) {
    if ( fabs(m_ChannelLowEnergy[i]-rmfChannelLowEnergy[i])/m_ChannelLowEnergy[i] > FUZZY ) return InconsistentChannels;
    if ( fabs(m_ChannelHighEnergy[i]-rmfChannelHighEnergy[i])/m_ChannelHighEnergy[i] > FUZZY ) return InconsistentChannels;
  }

  if ( m_EnergyUnits != r.getEnergyUnits() ) return InconsistentUnits;
  if ( m_RMFUnits != r.getRMFUnits() ) return InconsistentUnits;

  return OK;
}

Integer rmf::checkCompatibility(const arf& a) const
{
  vector<Real> arfLowEnergy = a.getLowEnergy();
  vector<Real> arfHighEnergy = a.getHighEnergy();
  if ( m_LowEnergy.size() != arfLowEnergy.size() ) return InconsistentEnergies;
  for ( size_t i=0; i < (size_t)m_LowEnergy.size(); i++ ) {
    if ( fabs(m_LowEnergy[i]-arfLowEnergy[i])/m_LowEnergy[i] > FUZZY ) return InconsistentEnergies;
    if ( fabs(m_HighEnergy[i]-arfHighEnergy[i])/m_HighEnergy[i] > FUZZY ) return InconsistentEnergies;
  }

  return OK;

}

// convert the energy units to keV if they are something else

Integer rmf::convertUnits()
{

  // set up energy/wave conversion factors and check for valid units

  bool xwave;
  Real xfactor;

  Integer status(OK);

  status = calcXfactor(m_EnergyUnits, xwave, xfactor);
  if ( status != OK ) return(status);

  if ( xfactor == 1.0 ) return(OK);

  if ( xwave ) {
    for (size_t i=0; i<m_LowEnergy.size(); i++) {
      Real temp(m_HighEnergy[i]);
      m_HighEnergy[i] = xfactor/m_LowEnergy[i];
      m_LowEnergy[i] = xfactor/temp;
    }
    for (size_t i=0; i<m_ChannelLowEnergy.size(); i++) {
      Real temp(m_ChannelHighEnergy[i]);
      m_ChannelHighEnergy[i] = xfactor/m_ChannelLowEnergy[i];
      m_ChannelLowEnergy[i] = xfactor/temp;
    }
  } else {
    for (size_t i=0; i<m_LowEnergy.size(); i++) {
      m_LowEnergy[i] = xfactor*m_LowEnergy[i];
      m_HighEnergy[i] = xfactor*m_HighEnergy[i];
    }
    for (size_t i=0; i<m_ChannelLowEnergy.size(); i++) {
      m_ChannelLowEnergy[i] = xfactor*m_ChannelLowEnergy[i];
      m_ChannelHighEnergy[i] = xfactor*m_ChannelHighEnergy[i];
    }
  }

  // if necessary reverse the rows in the response

  if (m_LowEnergy[0] > m_HighEnergy[m_HighEnergy.size()-1]) {
    this->reverseRows();
  }

  m_EnergyUnits = "keV";
  return(OK);
}

// reverse the rows, required if they are not in increasing order of energy
// note that this does not reverse the channels.
// note also that this does not the OrderGroup array

void rmf::reverseRows()
{

  // easiest just to set up a work rmf. Need to set the threshold and first channel in work rmf.

  rmf work;
  work.setResponseThreshold(m_ResponseThreshold);
  work.setFirstChannel(m_FirstChannel);

  // loop through the response in reverse order extracting the response vector
  // and constructing new response

  size_t Nbins(m_LowEnergy.size());
  vector<Real> ResponseValues(NumberChannels());

  for (size_t i=0; i<Nbins; i++) {

    ResponseValues = this->RowValues(Nbins-i-1);

    work.addRow(ResponseValues, m_LowEnergy[Nbins-i-1], m_HighEnergy[Nbins-i-1]);

  }

  // Reset the arrays based on those in work

  setFirstChannelGroup(work.getFirstChannelGroup());
  setNumberChannelsGroup(work.getNumberChannelsGroup());
  setMatrix(work.getMatrix());
  setNumberGroups(work.getNumberGroups());
  setLowEnergy(work.getLowEnergy());
  setHighEnergy(work.getHighEnergy());

  return;

}


void rmf::addRow(const vector<Real>& Response, const Real eLow, const Real eHigh) 
{

  Integer NGroups(0);
  bool inGroup(false);

  // temporary vectors for the new row

  vector<Integer> fchan(0), nchan(0);
  vector<Real> matrix(0);

  for ( size_t j=0; j<Response.size(); j++ ) {

    if ( Response[j] > m_ResponseThreshold ) {

      // if not in a response group then start a new one

      if ( !inGroup ) {

	NGroups++;

	fchan.push_back(j+m_FirstChannel);
	nchan.push_back(1);
	matrix.push_back(Response[j]);

	inGroup = true;

	// otherwise add next response to this group

      } else {

	nchan[nchan.size()-1]++;
	matrix.push_back(Response[j]);

      }

      // if response below threshold then end group if it is open

    } else {

      if ( inGroup ) inGroup = false;

    }

    // end loop over channels

  }

  // add the temporary vectors and number of groups and response energies

  m_FirstChannelGroup.push_back(fchan);
  m_NumberChannelsGroup.push_back(nchan);
  m_Matrix.push_back(matrix);
  m_NumberGroups.push_back(NGroups);
  m_LowEnergy.push_back(eLow);
  m_HighEnergy.push_back(eHigh);

  return;

}

// version of addRow for multiple grating orders

void rmf::addRow(const vector<vector<Real> >& Response, const Real eLow, const Real eHigh,
		 const vector<Integer>& GratingOrder) 
{

  Integer NGroups(0);
  bool First(true);

  // temporary vectors for the new row

  vector<Integer> fchan(0), nchan(0), order(0);
  vector<Real> matrix(0);

  // loop over grating orders

  for ( size_t i=0; i<Response.size(); i++ ) {

    bool inGroup(false);

    for ( size_t j=0; j<Response[i].size(); j++ ) {

      if ( Response[i][j] > m_ResponseThreshold ) {

	// if not in a response group then start a new one

	if ( !inGroup ) {

	  NGroups++;

	  fchan.push_back(j+m_FirstChannel);
	  nchan.push_back(1);
	  matrix.push_back(Response[i][j]);
	  order.push_back(GratingOrder[i]);

	  inGroup = true;

	  if ( First ) {
	    fchan.push_back(nchan.size()-1);
	    First = false;
	  }

	  // otherwise add next response to this group

	} else {

	  nchan[nchan.size()-1]++;
	  matrix.push_back(Response[i][j]);

	}

	// if response below threshold then end group if it is open

      } else {

	if ( inGroup ) inGroup = false;

      }

      // end loop over channels

    }

    // end loop over orders

  }

  // add the temporary vectors and the number of groups and response energies

  m_FirstChannelGroup.push_back(fchan);
  m_NumberChannelsGroup.push_back(nchan);
  m_Matrix.push_back(matrix);
  m_OrderGroup.push_back(order);
  m_NumberGroups.push_back(NGroups);
  m_LowEnergy.push_back(eLow);
  m_HighEnergy.push_back(eHigh);

  return;

}

// a version of addRow which uses channel groups

void rmf::addRow(const vector<Integer>& fChan, const vector<Integer>& nChan, 
		 const vector<Real>& Response, const Real eLow, const Real eHigh) 
{

  Integer NGroups(fChan.size());

  // add the vectors and number of groups and response energies

  m_FirstChannelGroup.push_back(fChan);
  m_NumberChannelsGroup.push_back(nChan);
  m_Matrix.push_back(Response);
  m_NumberGroups.push_back(NGroups);
  m_LowEnergy.push_back(eLow);
  m_HighEnergy.push_back(eHigh);

  return;

}

void rmf::substituteRow(const Integer RowNumber, const vector<Real>& Response) 
{

  // if RowNumber is invalid then don't do anything

  if ( RowNumber < 0 || RowNumber >= (Integer)m_NumberGroups.size() ) return;

  // construct compressed format arrays for the input response vector

  Integer NGroups(0);
  vector<Integer> FChan;
  vector<Integer> NChan;
  vector<Real> MatrixValues;
  bool inGroup(false);

  for ( size_t j=0; j<Response.size(); j++ ) {

    if ( Response[j] > m_ResponseThreshold ) {

      // if not in a response group then start a new one

      if ( !inGroup ) {

	NGroups++;

	FChan.push_back(j+m_FirstChannel);
	NChan.push_back(1);
	MatrixValues.push_back(Response[j]);

	inGroup = true;

	// otherwise add next response to this group

      } else {

	NChan[NChan.size()-1]++;
	MatrixValues.push_back(Response[j]);

      }

      // if response below threshold then end group if it is open

    } else {

      if ( inGroup ) inGroup = false;

    }

    // end loop over channels

  }

  // now replace the entries for the targeted row number

  m_FirstChannelGroup.at(RowNumber) = FChan;
  m_NumberChannelsGroup.at(RowNumber) = NChan;
  m_Matrix.at(RowNumber) = MatrixValues;

  // reset the number of groups for this row

  m_NumberGroups[RowNumber] = NGroups;

  return;

}

// version of substituteRow for multiple grating orders

void rmf::substituteRow(const Integer RowNumber, const vector<vector<Real> >& Response,
		 const vector<Integer>& GratingOrder) 
{

  // if RowNumber is invalid then don't do anything

  if ( RowNumber < 0 || RowNumber >= (Integer)m_NumberGroups.size() ) return;

  // construct compressed format arrays for the input response vector

  Integer NGroups(0);
  vector<Integer> FChan;
  vector<Integer> NChan;
  vector<Integer>GOrder;
  vector<Real> MatrixValues;

  // loop over grating orders

  for ( size_t i=0; i<Response.size(); i++ ) {

    bool inGroup(false);

    for ( size_t j=0; j<Response[i].size(); j++ ) {

      if ( Response[i][j] > m_ResponseThreshold ) {

	// if not in a response group then start a new one

	if ( !inGroup ) {

	  NGroups++;

	  FChan.push_back(j+m_FirstChannel);
	  NChan.push_back(1);
	  MatrixValues.push_back(Response[i][j]);
	  GOrder.push_back(GratingOrder[j]);

	  inGroup = true;

	  // otherwise add next response to this group

	} else {

	  NChan[NChan.size()-1]++;
	  MatrixValues.push_back(Response[i][j]);

	}

	// if response below threshold then end group if it is open

      } else {

	if ( inGroup ) inGroup = false;

      }
      
      // end loop over channels

    }

    // end loop over orders

  }

  // now replace the entries for the targeted row number

  m_FirstChannelGroup.at(RowNumber) = FChan;
  m_NumberChannelsGroup.at(RowNumber) = NChan;
  m_Matrix.at(RowNumber) = MatrixValues;
  m_OrderGroup.at(RowNumber) = GOrder;

  // reset the number of groups for this row

  m_NumberGroups[RowNumber] = NGroups;

  return;

}

// substitute a row into the response using channel group format

void rmf::substituteRow(const Integer RowNumber, const vector<Integer>& fChan,
			const vector<Integer>& nChan, const vector<Real>& Response) 
{

  // if RowNumber is invalid then don't do anything

  if ( RowNumber < 0 || RowNumber >= (Integer)m_NumberGroups.size() ) return;

  // replace the entries for the targeted row number

  m_FirstChannelGroup.at(RowNumber) = fChan;
  m_NumberChannelsGroup.at(RowNumber) = nChan;
  m_Matrix.at(RowNumber) = Response;

  // reset the number of groups for this row

  m_NumberGroups[RowNumber] = fChan.size();

  return;

}

// multiply a response by a vector and output a vector of pha values. The input
// vector is assumed to be on the energy binning

vector<Real> rmf::multiplyByModel(const vector<Real>& model)
{
  vector<Real> outPhaValues(0.0,this->NumberChannels());

  // loop over energies
  size_t nE = (size_t)(this->NumberEnergyBins());
  for (size_t ie=0; ie<nE; ie++) {

    Real modelValue(model[ie]);

    // loop over response groups for this energy

    size_t ir(0);
    for (size_t ig=0; ig<(size_t)m_NumberGroups[ie]; ig++) {

      for (size_t ich=(size_t)m_FirstChannelGroup[ie][ig]; 
	   ich<(size_t)(m_FirstChannelGroup[ie][ig]+m_NumberChannelsGroup[ie][ig]-1); ich++) {
	outPhaValues[ich] += modelValue * m_Matrix[ie][ir];
	ir++;
      }
    }
  }

  return outPhaValues;
}

// return a vector containing the FWHM in channels for each energy. This does
// assume that the response has a well-defined main peak and operates by the
// simple method of stepping out from the peak in both directions till the 
// response falls below half the maximum. A better solution would obviously be
// to fit a gaussian

vector<Real> rmf::estimatedFWHM() const
{
  size_t nE = (size_t)(this->NumberEnergyBins());

  vector<Real> fwhm(nE);

  for (size_t ie=0; ie<nE; ie++) {

    vector<Real> values = this->RowValues(ie);

    // find peak value
    Real maxValue = values[0];
    size_t imax = 0;
    for (size_t ich=1; ich<values.size(); ich++) {
      if ( values[ich] > maxValue ) {
	maxValue = values[ich];
	imax = ich;
      }
    }

    // now find the fwhm by moving outward from the maximum +ve and -ve directions
    // till we find the half maximum points or run into the edge.
    Real halfMax(maxValue/2.0);
    size_t ihigh = imax;
    if ( imax < values.size()-1 ) ihigh++;
    while ( ihigh < values.size()-1 && values[ihigh] > halfMax ) ihigh++;
    size_t ilow = imax;
    if ( ilow > 0 ) ilow--;
    while ( ilow > 0 && values[ilow] > halfMax ) ilow--;

    bool goodLow(true), goodHigh(true);
    if ( values[values.size()-1] > halfMax ) goodHigh = false;
    if ( values[0] > halfMax ) goodLow = false;

    fwhm[ie] = 0.0;
    if ( goodHigh) {
      fwhm[ie] += ihigh - imax;
    }
    if ( goodLow ) {
      fwhm[ie] += imax - ilow;
    }
    if ( (goodHigh && !goodLow) || (!goodHigh && goodLow) ) fwhm[ie] *= 2;
    if ( !goodHigh && !goodLow ) fwhm[ie] = -1.0;

  }

  return fwhm;
}

// return a vector containing the FWHM in channels for each channel. This does
// assume that the response has a well-defined main peak

vector<Real> rmf::estimatedFWHMperChannel() const
{
  size_t nE = (size_t)(this->NumberEnergyBins());
  size_t nChan = (size_t)(this->NumberChannels());

  vector<Real> Efwhm(nE);
  vector<Real> fwhm(nChan);

  // first estimate the FWHM for each energy bin

  Efwhm = this->estimatedFWHM();

  // now interpolate using the nominal channel energies to give the FWHM for
  // each channel. assuming that FWHM does not change significantly over the
  // channel so just find the FWHM at the center energy of the channel

  for (size_t i=0; i<nChan; i++) {

    Real channelE = 0.5*(m_ChannelLowEnergy[i] + m_ChannelHighEnergy[i]);
    size_t index = binarySearch(channelE, m_LowEnergy, m_HighEnergy);

    fwhm[i] = Efwhm[index];

  }

  return fwhm;
}


// utility routines that are not methods for the rmf object

rmf operator* (const rmf& r, const arf& a){
  rmf rr(r);
  rr *= a;
  return rr;
}

rmf operator* (const arf& a, const rmf& r){
  rmf rr(r);
  rr *= a;
  return rr;
}

rmf operator* (const rmf& r, const Real& f){
  rmf rr(r);
  rr *= f;
  return rr;
}

rmf operator* (const Real& f, const rmf& r){
  rmf rr(r);
  rr *= f;
  return rr;
}

rmf operator+ (const rmf& a, const rmf& b){
  rmf r(a);
  r += b;
  return r;
}

// return an integer vector containing the extension numbers of all EBOUNDS
// extensions in the file

vector<Integer> RMFeboundsExtensions(const string filename, const string extname)
{
  vector<Integer> extNums(0);

  // Attempt to open filename as a FITS object

  unique_ptr<FITS> pInfile((FITS*)0);
  try {
    pInfile.reset(new FITS(filename,Read));
  } catch(...) {
    string msg = "Failed to open "+filename;
    SPreportError(NoSuchFile, msg);
    return extNums;
  }

  // loop round the extensions finding the EBOUNDS or the value input as
  // extname

  string findName("EBOUNDS");
  if ( extname.length() != 0 ) findName = extname;

  int i = 1;
  while (true) {
    try {
      ExtHDU& test = pInfile->extension(i);
      string defValue = "UNKNOWN";
      if ( SPreadKey(test,"EXTNAME",defValue) == findName ) {
	extNums.push_back(i);
      } else {
	if ( SPreadKey(test,"HDUCLAS1",defValue) == "RESPONSE" &&
	     SPreadKey(test,"HDUCLAS2",defValue) == "EBOUNDS" ) {
	  extNums.push_back(i);
	}
      }
      i++;
    } catch(...) {
      return extNums;
    }
  }
  return extNums;
}

// return an integer vector containing the extension numbers of all MATRIX
// extensions in the file

vector<Integer> RMFmatrixExtensions(const string filename, const string extname)
{
  vector<Integer> extNums(0);

  // Attempt to open filename as a FITS object

  unique_ptr<FITS> pInfile((FITS*)0);
  try {
    pInfile.reset(new FITS(filename,Read));
  } catch(...) {
    string msg = "Failed to open "+filename;
    SPreportError(NoSuchFile, msg);
    return extNums;
  }

  // loop round the extensions finding the MATRIX or if extname input is 
  // non-zero length use that

  string findName("MATRIX");
  if ( extname.length() != 0 ) findName = extname;

  int i = 1;
  while (true) {
    try {
      ExtHDU& test = pInfile->extension(i);
      string defValue = "UNKNOWN";
      if ( SPreadKey(test,"EXTNAME",defValue) == findName ) {
	extNums.push_back(i);
      }	else if ( SPreadKey(test,"EXTNAME",defValue) == "SPECRESP MATRIX" ) {
	extNums.push_back(i);
      } else {
	if ( SPreadKey(test,"HDUCLAS1",defValue) == "RESPONSE" &&
	     SPreadKey(test,"HDUCLAS2",defValue) == "RSP_MATRIX" ) {
	  extNums.push_back(i);
	}
      }
      i++;
    } catch(...) {
      return extNums;
    }
  }
}

// calculate the response vector for some energy given a gaussian width
// the gaussian is assumed to be in the units of energy,
// ChannelLowEnergy and ChannelHighEnergy

void calcGaussResp(const Real sigma, const Real energy, const Real threshold, 
		   const vector<Real>& channelLowEnergy, 
		   const vector<Real>& channelHighEnergy, 
		   vector<Real>& ResponseVector)
{

  Real winv = 1.0/sigma/sqrt(2.0);
  size_t N = channelLowEnergy.size();
  ResponseVector.resize(N);
  for (size_t i=0; i<N; i++) ResponseVector[i] = 0.0;

  // find the channel containing the energy

  size_t icen = binarySearch(energy, channelLowEnergy, channelHighEnergy);

  // first do the case of zero line width

  if ( sigma <= 0.0 ) {
    ResponseVector[icen] = 1.0;
    return;
  }

  // if the line center is below the first bin then don't calculate the lower
  // part of the line. If the line center is above the last bin then just calculate
  // the part of the line within the energy range

  if ( energy < channelLowEnergy[0] ) {
    icen = 0;
  } else if ( energy > channelHighEnergy[N-1] ) {
    icen = N-1;
  }

  // Do the low energy part of the line

  Integer ielow((Integer)icen);
  Real alow(0.0);
  Real lineSum(0.0);
  Real ahi;

  while ( ielow >= 0 ) {
    ahi = erf(winv*(fabs(channelLowEnergy[ielow]-energy)));
    Real fract = (ahi-alow)/2;
    if ( fract >= threshold || ielow == (Integer)icen ) {
      ResponseVector[ielow] += fract;
      lineSum += fract;
    } else {
      ielow = 0;
    }
    alow = ahi;
    ielow -= 1;
  }

  // If the line center is above the last bin then don't calculate the upper
  // part of the line. If line center is below the first bin then just calculate
  // the part of the line within energy range

  if ( energy < channelLowEnergy[0] ) {
    icen = 1;
  } else if ( energy > channelHighEnergy[N-1] ) {
    icen = N + 1;
  }

  // Do the high energy part of the line

  ielow = icen;
  alow = 0.0;
  while ( ielow <= (Integer)N-1 ) {
    ahi = erf(winv*(fabs(channelHighEnergy[ielow]-energy)));
    Real fract = (ahi-alow)/2;
    if ( fract >= threshold || ielow == (Integer)icen ) {
      ResponseVector[ielow] += fract;
      lineSum += fract;
    } else {
      ielow = N;
    }
    alow = ahi;
    ielow += 1;
  }

  return;
}

size_t binarySearch(const Real energy, const vector<Real>& lowEnergy,
		    const vector<Real>& highEnergy)
{
  // Function to do a binary search for the i which satisfies
  // lowEnergy[i] < energy <= highEnergy[i]
  
  size_t nE = lowEnergy.size();
  
  bool increase(false);
  if ( lowEnergy[1] > lowEnergy[0] ) increase = true;
  
  if ( increase ) {
    if ( energy < lowEnergy[0] ) return(0);
    if ( energy > highEnergy[nE-1] ) return(nE-1);
  } else {
    if ( energy > lowEnergy[0] ) return(0);
    if ( energy < highEnergy[nE-1] ) return(nE-1);
  }
  
  size_t low = 0;
  size_t high = nE-1;
  size_t bisearch;
  
  while ( high-low > 1 ) {
    bisearch = (low+high)/2;
    if ( (increase && energy > lowEnergy[bisearch]) ||
	 (!increase && energy < lowEnergy[bisearch]) ) {
      low = bisearch;
    } else {
      high = bisearch;
    }
  }

  if ( lowEnergy[low] < energy && energy <= highEnergy[low] ) {
    return(low);
  } else {
    return(high);
  }

}