//   Read the documentation to learn more about C++ code generator
//   versioning.
//	  %X% %Q% %Z% %W%

// DataContainer
#include <XSModel/GlobalContainer/DataContainer.h>
// Grid
#include <XSFit/Fit/Grid.h>
// SpectralData
#include <XSModel/Data/SpectralData.h>
// Graph
#include <XSPlot/Plot/Graph.h>
// Background
#include <XSModel/Data/BackCorr/Background.h>
// ModelContainer
#include <XSModel/GlobalContainer/ModelContainer.h>
// Model
#include <XSModel/Model/Model.h>
// Fit
#include <XSFit/Fit/Fit.h>
// Plot
#include <XSPlot/Plot/Plot.h>

#include <XSFit/StatMethod/ChiSquare/ChiSquare.h>
#include <XSPlot/Plot/GraphCreator.h>
#include <XSPlot/Plot/PlotCommandCreator.h>
#include <XSPlot/Plot/PlotGroup.h>
#include <XSModel/Model/Component/SumComponent.h>
#include <XSUtil/Parse/XSparse.h>
#include <XSUtil/Numerics/Numerics.h>
#include <XSsymbol.h>
#include <XSstreams.h>
#include <XSContainer.h>
#include <functional>
#include <fstream>
#include <sstream>
#include <iterator>
#include <set>
#include <cmath>
#include <algorithm>

Plot*   XSContainer::plot = 0;
using namespace XSContainer;


// Class Plot::NoEnergiesDefined 

Plot::NoEnergiesDefined::NoEnergiesDefined (const string& diag)
   : YellowAlert(" no energies defined for spectrum ")
{
  tcerr << diag << '\n';
}


// Class Plot::NoDataLoaded 

Plot::NoDataLoaded::NoDataLoaded()
  : YellowAlert(" no data loaded\n")
{
}


// Class Plot::NoModelDefined 

Plot::NoModelDefined::NoModelDefined (const string& command)
  : YellowAlert(" no model defined\n")
{
  tcerr << command << '\n';
}


// Class Plot::ModelNotFound 

Plot::ModelNotFound::ModelNotFound (const string& msg)
  : YellowAlert(" model not found\n")
{
  tcerr << msg << '\n';
}


// Class Plot::RebinInfo 

// Class Plot::NoGrid 

Plot::NoGrid::NoGrid (const string& diag)
  : YellowAlert(" no grid defined: ")
{
  tcerr << diag << '\n';
}


// Class Plot::InvalidContourLevels 

Plot::InvalidContourLevels::InvalidContourLevels (const string& diag)
  : YellowAlert(" contour level specification is invalid: ")
{
  tcerr << diag << '\n';
}


// Class Plot::NotImplemented 

Plot::NotImplemented::NotImplemented (const string& diag)
  : YellowAlert(" plot type not implemented: ")
{
  tcerr << diag << '\n';
}


// Class Plot 
Plot* Plot::s_instance = 0;
const Real Plot::s_NODATA = -1.2e-34;
const int Plot::s_MAXCMDLEN = 72;

Plot::Plot (const string& graphPackage, const string& labelDictionaryPath)
      : m_isCurrent(false),
        m_showAddComponents(false),
	m_xOption(CHANNELS),
	m_prevGroupRange(0,0),
	m_plotIDs(false),
	m_temperature(1.0),
	m_emisLimit(1.0e-18),
	m_redshift(0.),
	m_labelDictionaryFile("plotLabels.dat"),
        m_counts(false),
	m_selectedLabels(3,""),
        m_yLog(false),
	m_rangeDefaultFlags(0),
	m_xLog(true),
	m_dataStyle(PlotStyle::DOT),
	m_modelStyle(PlotStyle::SOLID),
	m_xLine(false,.0),
	m_plotXLog(false),
	m_plotYLog(false),
	m_numGraphs(1),
	m_currentGraphNumber(0),
	m_modelStatus(INACTIVE),
	m_plotBeginIndex(0),
	m_groupsRebinInfo(),
	m_lastRebinEntry(),
        m_plotType(VECTOR),
	m_lastCommand(1, "data"),
        m_grid(0),  // non-owning pointer
        m_divideByArea(false,false),
        m_savedPlotArray(),
        m_saveArrayInfo(Plot::NO_SAVE, 1),
        m_splashPage(true),
        m_modelPositionDirectory(),
        m_modelsForSpectra(),
        m_showBackground(false),
        m_waveUnitsHz(1),
        m_ranges(),
        m_spectra(),
        m_labelNames(),
        m_graph(0),
        m_groups()
{
       // Register the Plot object as an observer of the DataContainer.
       // This will create the DataContainer if it doesn't already exist -
       // safety measure.
       XSContainer::datasets = XSContainer::DataContainer::Instance();
       XSContainer::datasets->Attach(this);

       m_ranges["x"] = std::pair<Real,Real>(0.,1.);
       m_ranges["y"] = std::pair<Real,Real>(0.,1.);
       m_ranges["z"] = std::pair<Real,Real>(0.,1.);

       RebinInfo initEntry;
       initEntry.sigma = 0; initEntry.maxBins = 1; initEntry.mode = STD;
       m_lastRebinEntry = std::pair<int, RebinInfo>(-1, initEntry);
       m_groupsRebinInfo.insert(m_lastRebinEntry);

  loadLabelDictionary(labelDictionaryPath);
  m_graph = GraphCreator::MakeGraph(graphPackage);
}


Plot::~Plot()
{
   PlotGroupContainer::iterator pG = m_groups.begin();
   while ( pG != m_groups.end())
   {
        delete *pG; 
        ++pG;      
   }
   m_groups.clear();
   delete m_graph;
   PlotCommandCreator::destroy();
}


Plot* Plot::Instance (const string& graphPackage, const string& labelDictionaryPath)
{
  if (s_instance == 0) 
  {
     s_instance = new Plot(graphPackage, labelDictionaryPath);
  } 
  return s_instance;
}

void Plot::Update (Subject* changed)
{
  initializeSpectra(static_cast<XSContainer::DataContainer*>(changed));
}

void Plot::setXOption (const string& option)
{
  if (option == "channel")
  {
     m_xOption = CHANNELS;
     SpectralData::energySpace(false);
     SpectralData::wavelengthSpace(false);    
  }
  else if (option == "energy")
  {
     m_xOption = ENERGY;
     SpectralData::energySpace(true);
     SpectralData::wavelengthSpace(false);
  }
  else if (option == "wavelength")
  {
     m_xOption = WAVELENGTH;
     SpectralData::energySpace(false);
     SpectralData::wavelengthSpace(true);
  }
  else throw RedAlert(" Programmer Error: X-axis mode");
}

void Plot::makeBinnedPlotArrays (bool unfold)
{
  size_t Ngroups (datasets->numberOfPlotGroups()); 
  itemizePlotGroupModels();

  std::vector<std::list<SpectralData*> > groupedSpectra(Ngroups,std::list<SpectralData*>());

  m_ranges["x"].first =  1.E+100;    
  m_ranges["x"].second = -1;    

  for (size_t j = 0; j <= Ngroups - 1; ++j )
  {
     int chanXGroupOffset (0); // icbeg
     // gtdpin...    
     PlotGroup* gr = 0; 
     // can throw an "EnergiesNotDefined" exception...
     gr = initializePlotGroup(j+1,groupedSpectra[j],unfold);
     if (j == 0) m_plotBeginIndex = m_groups.size() - 1;
     if (gr && gr->n > 0)
     {
        int npts (0);
        int latestChannel(0); // ice
        bool moreBins (true);
        try 
        {
           if (m_divideByArea.second) initEffectiveAreas(groupedSpectra[j]);
           do // execute at least once...
           {
              if ( moreBins = getNextPlotBin(gr, npts, groupedSpectra[j], 
                                      latestChannel,unfold))
              {
                 if ( m_xOption == CHANNELS ) gr->xAxis.data[npts] += chanXGroupOffset;
                 ++npts;
              }
	      else gr->n = npts;

           }  while ( moreBins );
           if (m_divideByArea.second) clearEffectiveAreas(groupedSpectra[j]);
        }
        catch (...)
        {
           // Whether or not initEffectiveAreas was called,
           // this is safe to do:
           clearEffectiveAreas(groupedSpectra[j]);
           resetPlotGroups();
           throw;
        }
        if (gr->model.size() > 1)
        {
	   std::vector<Real>& modelSum = gr->model[0].data;
           for (size_t l = 1; l < gr->model.size(); ++l)
           {
	      std::vector<Real>& model = gr->model[l].data;
              for (int k = 0; k < npts; ++k) 
              {
                      modelSum[k] += model[k];       
              }
           }
        }
        Real xAxis_0 = gr->xAxis.data[0];
        Real xAxis_N = gr->xAxis.data[npts-1];
        if (xAxis_N > xAxis_0)
        {
           m_ranges["x"].first = std::min(xAxis_0
                        - gr->xAxis.errors[0][0],m_ranges["x"].first);
           m_ranges["x"].second  = std::max(xAxis_N 
                        + gr->xAxis.errors[0][npts-1],m_ranges["x"].second);
        }
        else
        {
           m_ranges["x"].first = std::min(xAxis_N 
                        - gr->xAxis.errors[0][npts-1],m_ranges["x"].first);                                
           m_ranges["x"].second = std::max(xAxis_0 
                        + gr->xAxis.errors[0][0],m_ranges["x"].second);
        }
        chanXGroupOffset += gr->channelPlotSpacing;
     }
  }
}

void Plot::makePlotArraysFromSpectra ()
{
  // makes one PlotGroup object per spectrum response
  // and ignore the plot grouping mechanism.

  size_t NS (datasets->numberOfSpectra());

  if ( m_xOption == CHANNELS || m_xOption == ENERGY )
  {
     m_ranges["x"].first =  1.E+100;    
     m_ranges["x"].second = -1;    
  }
  else
  {
     m_ranges["x"].first = -1;
     m_ranges["x"].second = 1.E+100;             
  }  
  for ( size_t s = 1; s <= NS; ++s )
  {
     SpectralData* sp(datasets->lookup(s));
     try
     {
        initializePlotGroup(sp);
     }
     catch (Plot::NoEnergiesDefined)
     {
     }        
  } 

  if (m_groups.empty())
     throw YellowAlert("No response energies for any loaded spectra.\n");
}

void Plot::setPlotGroupNums (const string& rangesString)
{

    size_t sz (datasets->numberOfSpectra());
    if (!sz)
    {
       tcout << "No spectra loaded, cannot assign plot group numbers."<<std::endl;
       return;
    }
    if (!m_prevGroupRange.first)  m_prevGroupRange.first = 1;
    if (!m_prevGroupRange.second) m_prevGroupRange.second = sz;
    std::pair<size_t,size_t> limits = std::make_pair<size_t,size_t>(1,sz);

    string::size_type iSpace = rangesString.find_first_of(' ',0);
    string::size_type prevSpace = 0;
    while(prevSpace != string::npos)
    {
	string subString = rangesString.substr(prevSpace, iSpace-prevSpace);
	std::pair<size_t,size_t> range = 
                XSparse::wildRange(subString, limits, m_prevGroupRange);
        size_t begin = range.first;
        size_t end = range.second;

	if (begin > sz || end > sz)
	{
	    std::ostringstream msg;
	    if (begin == end)
	    {
		msg <<'\n'<< begin <<" is out of range of existing spectra";
	    }
	    else
	    {
		msg <<'\n' << begin << '-' << end 
		    << " extends out of range of existing spectra";
	    }   
	    throw XSparse::InvalidRange(msg.str());
	}

	// Begin and end are 1-based.

	// Check that the bins match up before attempting to group
	// spectra together.
	bool isOKtoGroup = true;
	SpectralData* sd=datasets->lookup(begin);

	for (size_t i=begin + 1; i <= end; ++i)
	{
	    SpectralData* sdata = datasets->lookup(i);

	    if (!sdata->energiesEqual(sd))
	    {
		tcout << "***Warning:  Spectra " << begin << '-' << end
		      << " do not contain the same energy bins and/or channels.\n"
		      << "   They will not be grouped.\n";
		isOKtoGroup = false; 
		break;
	    }
	}

	if (isOKtoGroup)
	{     
	    // Find highest group number prior to spectrum[begin-1].
	    int high = 0;
	    for (size_t i = 0; i < begin-1; ++i)
	    {
		if (datasets->plotGroupNums(i) > high)
		{
		    high = datasets->plotGroupNums(i);
		}
	    }     
	    // Assign spectra begin-1 to end-1 their new group number.
	    ++high;
	    for (size_t i = begin; i <= end; ++i)
	    {
		datasets->setPlotGroupNums(i-1, high);
	    }     
	    // Now correct the later spectra accordingly.
	    datasets->renumberPlotGroups(end+1, high);
	}     
	// prevSpace is actually set here to the next char after the space.
	prevSpace = rangesString.find_first_not_of(' ', iSpace);
	iSpace = rangesString.find_first_of(' ', prevSpace);
    }
    initializeSpectra(datasets);
}

void Plot::setIDs (const std::vector<Real>& values)
{

   m_plotIDs = true;
   size_t sz = values.size();
   if (sz)  m_temperature = values[0];
   if (sz > 1)  m_emisLimit = values[1];
   if (sz > 2)  m_redshift = values[2];
}

void Plot::loadLabelDictionary (const string& dictionaryPath)
{


   string fullPathName = dictionaryPath + string("/") + m_labelDictionaryFile;
   std::ifstream labelFile;
   labelFile.open(fullPathName.c_str());

   if (!labelFile)
   {
      string message = "Cannot open label dictionary file: " + 
      			m_labelDictionaryFile;
      throw RedAlert(message);
   }

   while (!labelFile.eof())
   {
      string line;
      getline(labelFile, line);
      size_t sz = line.size();
      if (sz != 0)
      {
         size_t i;
         string key("");
	 // key will equal the first sequential group of non-whitespace
	 // characters encountered in the line.
	 for (i=0; i<sz; ++i)
	 {
	    if (!isspace(line[i]))  break;
	 }
	 while (i < sz)
	 {
	    if (isspace(line[i]))  break;
	    key += line[i];
	    ++i;
	 }

	 if (key.size() != 0)
	 {
	    // Now proceed to the next non-whitespace character then read from
	    // it until the end of the line into string val.
	    string val("");
	    while (i < sz)
	    {
	       if (!isspace(line[i]))
	       {
	          val = line.substr(i,sz-i);
		  break;
	       }
	       ++i;
	    }
	    // Remove trailing whitespace from end of val and insert into map.
	    int j;
	    for (j=val.size()-1; j>=0; --j)
	    {
	       if (!isspace(val[j]))  break;
	    }
	    if (j >= 0)
	    {
	       m_labelNames.insert(std::pair<string,string>
			  (key, string(val,0,j+1))); 
	    }
	 }

      }
   }
   if (m_labelNames.size() == 0)
   {
      string message = "No valid entries found in " + m_labelDictionaryFile;
      throw RedAlert(message);
   }
}

bool Plot::getNextPlotBin (PlotGroup* gr, int npts, std::list<SpectralData*>& groupSpectra, int& latestChannel, bool unfold)
{

  if (groupSpectra.empty() || latestChannel > gr->lastInputChannel  ) return false;

  // fortran has 3 integer output variables, icb, ice, nbinrt.
  // nbinrt is only used as a boolean value and is returned by this function.
  // icb is the bin number in the input spectra array of the beginning of
  // the current bin, and ice the end of the current bin.

  Real chminl(0), chmaxl(0), chminh(0), chmaxh(0);
  Real count (0), backgroundCount(0);

  // if there is more than one source per spectrum, the models are stored
  // in the model std:vector<Real*> container enumerated from 1.
  // if there isn't then that vector only has one element, enumerated as zero.
  // for the latter case, the source number is 1 and the binned model is stored
  // in location zero. for the former case, the binned models are stores in
  // location == source number and the sum will ultimately be stored in location
  // zero.

  bool doneBin (false);
  bool notChannelSpace(m_xOption != CHANNELS);

  int chan (latestChannel);
  int firstChannelInBin = -1;
  int chanEnd (gr->lastInputChannel);
  int nBeg (0); // number of input bins making up the first output bin
  int numChannelsBinned(0); 
  Real critVar ( gr->critSignificance * gr->critSignificance );  
  Real& datapt = gr->yData.data[npts];
  Real& errorpt = gr->yData.errors[0][npts];
  Real& backpt = gr->background.data[npts];
  Real& backErrorpt = gr->background.errors[0][npts];
  // Using saveData array as a miscellaneous storage area.  By default
  // it will pass along area*time info for the bin, and is needed
  // for plot chisq when using C-stat.
  Real& areaSum = gr->saveData[npts];
  Real backSum = .0;
  Real backScaleRatio = .0;
  Real nBackChans = 0;
  const std::map<size_t,size_t>& modAssignments =
        m_modelPositionDirectory[gr->plotGroup-1];

  while (!doneBin )
  {
     int nEnd (0);
     Real tchmaxl(0);
     Real tchmaxh(0);

     std::list<SpectralData*>::const_iterator gs    (groupSpectra.begin());
     std::list<SpectralData*>::const_iterator gsEnd (groupSpectra.end());

     // gr->bin is the point in the SpectralData array under consideration.
     // chan is the current cha
     // npts is the current bin of the output array.
     while (gs != gsEnd && chan < static_cast<int>(gr->n))
     {       
        const SpectralData* current = *gs;
        size_t spectrumNumber (current->spectrumNumber());
        int offset = current->startChan() - current->firstChan();
        // if we got this far, detector must exist or be irrelevant.
        const Response* firstResp = 0;
        std::vector<Response*>::const_iterator itResp = current->detector().begin();
        while (itResp != current->detector().end() && !firstResp)
        {
           if (*itResp)
              firstResp = *itResp;
           ++itResp;
        }
        if ( (*gs)->noticedChannels(chan + offset) )
        {
           if ( firstChannelInBin == -1 )
           {
              firstChannelInBin = chan;
              chanEnd = std::min(chanEnd,gr->critNumChans+chan-1);
           }

           if ( chan == firstChannelInBin )
           {
              ++nBeg;
              if ( notChannelSpace )
              {
                 chminl += firstResp->eboundsMin()[chan];
                 chminh += firstResp->eboundsMax()[chan];
              }
           }

           if ( notChannelSpace )
           {
              tchmaxl += firstResp->eboundsMin()[chan];
              tchmaxh += firstResp->eboundsMax()[chan];
           }

           ++nEnd;
           ++numChannelsBinned;

           // accumulate data from SpectralData
           Real areaTimeFactor ( current->exposureTime()*current->areaScale(chan));
           areaSum += areaTimeFactor;
           Real factor = m_counts ? areaTimeFactor : 1;
           if (m_divideByArea.second)
           {
              Real effArea = current->effectiveAreas()[chan];
              if (effArea == 0.0)
              {
                 std::ostringstream msg;
                 msg << "Cannot divide data by area.\n   Zero effective area found in spectrum "
                    << spectrumNumber << " channel " << chan;
                 throw YellowAlert(msg.str());
              }
              factor /= effArea;
           }
           datapt += current->spectrum(chan)*factor; 
           count += current->spectrum(chan)*areaTimeFactor;
           errorpt += factor*factor*current->variance(chan);
           if ( current->background() )
           {
              const SpectralData* bckSpec = current->background()->data();
              const RealArray& backgroundSpectrum = bckSpec->spectrum();
              const RealArray& backgroundVariance = bckSpec->variance();
              Real backContribution = backgroundSpectrum[chan]*factor;
              backpt += backContribution;
              datapt -=  backContribution;
              backgroundCount += backgroundSpectrum[chan]*areaTimeFactor;
              Real backErrContribution = factor*factor*backgroundVariance[chan];
              errorpt += backErrContribution;
              backErrorpt += backErrContribution;
              // The following info is only needed if a minimum
              // variance needs to be applied in setError.  Taking
              // the average of these quantities across the bins is
              // merely a convenient way (and not necessarily the best)
              // for determining a minimum variance rate.
              backSum += bckSpec->exposureTime()*bckSpec->areaScale(chan);
              backScaleRatio += 
                current->backgroundScale(chan)/bckSpec->backgroundScale(chan);
              ++nBackChans;

           }
           if ( current->correction())
           {
              const RealArray& correctionSpectrum 
                              = current->correction()->spectrum();
              datapt -= current->correctionScale()*
                                      correctionSpectrum[chan]*factor;
              count -=  current->correctionScale()*
                              correctionSpectrum[chan]*areaTimeFactor;      
           }   


           std::vector<const Model*>::const_iterator itMod = 
                m_modelsForSpectra[spectrumNumber-1].begin();
           std::vector<const Model*>::const_iterator itModEnd = 
                m_modelsForSpectra[spectrumNumber-1].end();
           while (itMod != itModEnd)
           {
              const Model* m = *itMod;
              const size_t sourceIndex = 
                        modAssignments.find(m->sourceNumber())->second;
              const ArrayContainer& folded = m->foldedModel();
              ArrayContainer::const_iterator mmf = folded.find(spectrumNumber);
              if (mmf != folded.end())
              {
                  const RealArray& spectrumFit = mmf->second;
                  gr->model[sourceIndex].data[npts] 
                                    += spectrumFit[chan]*factor;                                           
              }
              if ( m_showAddComponents || unfold )
              {
                  // > 1: because we don't want to do this 
                  // if there is only one component source
                  const SourceList& sources = m->sources();
                  const size_t compOffset = gr->model.size() > 1 ? 1 : 0;
                  if (!gr->sources[sourceIndex-compOffset].empty())
                  {
                     SourceList::const_iterator sl = sources.begin();    
                     SourceList::const_iterator slEnd = sources.end(); 
                     PlotVectorList::iterator vl  = 
                        gr->sources[sourceIndex-compOffset].begin();
                     PlotVectorList::iterator vfl = 
                        gr->sourceFluxes[sourceIndex-compOffset].begin();
                     while (sl != slEnd)
                      {                          
                         ArrayContainer::const_iterator slf 
                                  = (*sl)->foldedPhotonFlux().find(spectrumNumber);
                         if ( slf != (*sl)->foldedPhotonFlux().end())
                         {
                            const RealArray& sourceFit = slf->second;
                            vl->data[npts] += sourceFit[chan]*factor;       
                         }
                         if (unfold)
                         {
                              ArrayContainer::const_iterator slff 
                                  = (*sl)->photonArray().find(spectrumNumber);
                              if ( slff != (*sl)->photonArray().end())
                              {
                                 const RealArray& sourceFit = slff->second;
                                 vfl->data[npts] += sourceFit[chan]*factor;       
                              }
                              ++vfl;
                         }
                         ++vl, ++sl;
                      }  
                  } 
              } // end if show/add components
              ++itMod;
           }
           latestChannel = chan + 1;
        } // end if noticed channel                   
        ++gs;
     } // end loop over spectra in group

     // Calculate chanEmax for the last bin with a non-zero
     // contribution to the accumlated sums.                
     if (nEnd > 0)
     {
        chmaxl = tchmaxl/nEnd;       
        chmaxh = tchmaxh/nEnd;       
     }

     doneBin = ( chan >= chanEnd ) || ( datapt > 0 && datapt*datapt >= errorpt*critVar);

     ++chan;                        
  } // end while !doneBin

  if ( numChannelsBinned == 0 ) return false;

// Set the channel min energies

      if ( nBeg > 0 ) 
      {
         chminl = chminl/nBeg;
         chminh = chminh/nBeg;
      }

// Now set the min and max channel energies. This convoluted procedure is
// necessary to handle the cases of both channels increasing and decreasing
// with energy

  Real chanEmin = std::min ( std::min(chminl,chminh), std::min(chmaxl, chmaxh) );
  Real chanEmax = std::max ( std::max(chminl,chminh), std::max(chmaxl, chmaxh) );

  areaSum /= numChannelsBinned;
  if (nBackChans)
  {
     backSum /= nBackChans;
     backScaleRatio /= nBackChans;
  }
  else
  {
     // Negative area*time tells calcMinVariance function that no
     // background exists.
     backSum = -1.0; 
     backScaleRatio = 1.0;
  }

// Set the error. If we are summing in quadrature simply convert from 
// variance to standard deviations.
// N.B. Should check the error method used against the variance weighting
// scheme set for data (DataUtility::statWeight). 

  setError( gr->errorType, errorpt,  backErrorpt, areaSum, backSum, count, backgroundCount, backScaleRatio );

// Set the scaling factor to get data, error, and model in the correct
// units (note that the energy and wavelength normalizations assume that
// the no. of square cm. is the same in each grouped file). 


  Real norm(0);

  switch ( m_xOption)
  {
     default:
     case CHANNELS:
         norm = 1./Real(numChannelsBinned);
         break;
     case ENERGY:
        {
           Real chans(latestChannel - firstChannelInBin);
           norm = (chanEmin == chanEmax) ? 1 
                 : chans/(std::abs(chanEmax-chanEmin)*Real(numChannelsBinned));            
        }    
        break;
     case WAVELENGTH:
        {
            Real chans(latestChannel - firstChannelInBin);
            if (m_waveUnitsHz == 1)
               norm =  ( chanEmin == chanEmax || chanEmin <= 0 || chanEmax <= 0 ) ? 1
                   : chans/(std::abs( Numerics::KEVTOHZ*(chanEmax - chanEmin)
                                             *Real(numChannelsBinned)));
            else
               norm =  ( chanEmin == chanEmax || chanEmin <= 0 || chanEmax <= 0 ) ? 1
                   : chans/(std::abs( Numerics::KEVTOA*(1./chanEmax - 1./chanEmin)
                                             *Real(numChannelsBinned)));

        }    
        break;    
   }

   if ( norm != 1)
   {
        errorpt *= norm;
        datapt  *= norm;  
        backErrorpt *= norm;
        backpt *= norm;
        if ( gr->model.size() == 1) 
        {
           gr->model[0].data[npts] *= norm;
        }
        else
        {
           for ( size_t k = 1;  k < gr->model.size(); ++k)     
           {
              gr->model[k].data[npts] *= norm;       
           }
        }

        if ( m_showAddComponents || unfold )
        {
            for (size_t j = 0 ; j  < gr->sources.size(); ++j )
            {   
                PlotVectorList::iterator vl = gr->sources[j].begin();    
                PlotVectorList::iterator vlEnd = gr->sources[j].end(); 
                while (vl != vlEnd)
                {                          
                   vl->data[npts] *= norm;
                   ++vl;
                }   
                if (unfold)
                {
                   PlotVectorList::iterator fl = gr->sourceFluxes[j].begin();    
                   PlotVectorList::iterator flEnd = gr->sourceFluxes[j].end(); 
                   while (fl != flEnd)
                   {                          
                      fl->data[npts] *= norm;
                      ++fl;
                   }   
                }
            }
        }        
    }

    // flux arrays

    if (unfold)
       integrateModelFluxes(groupSpectra,chanEmin,chanEmax,gr,npts);

    // x axis and x bin width arrays. could be done earlier but better
    // to separate responsibilities.

    setXaxis(gr,firstChannelInBin,latestChannel,chanEmin,chanEmax,npts);

  return ( numChannelsBinned > 0 );
}

PlotGroup* Plot::initializePlotGroup (size_t groupIndex, std::list<SpectralData*>& groupSpectra, bool unfold)
{

    // first: get the spectra which belong to this plot group.
    // all of the spectra in the group have been checked to have the same number of
    // channels.      


    std::pair<SpecGroupIter,SpecGroupIter> r (m_spectra.equal_range(groupIndex));     
    SpecGroupIter s ( r.first );
    bool notChannels = (m_xOption != CHANNELS);


    while ( s != r.second )
    {
        SpectralData* current (s->second);
        // test for "ignored" data and do not add to group if no channels are there.
        const BoolArray& noticed = current->noticedChannels();
        BoolArray::const_iterator t (std::find(noticed.begin(),noticed.end(),true));

        if ( t != noticed.end())
        {
           groupSpectra.push_back(current); 
           if (notChannels)
           {
              std::vector<Response*>::const_iterator itDet = current->detector().begin();         
              std::vector<Response*>::const_iterator itDetEnd = current->detector().end();
              while (itDet != itDetEnd)
              {
                 if (*itDet && (*itDet)->eboundsMin().size())
                    break;
                 ++itDet;
              }
              if (itDet == itDetEnd)         
              {
	          std::ostringstream diag;
	          diag << current->spectrumNumber() << ", but plot requested on "
		       << "energy/wavelength axis.\n Plot Group " << groupIndex 
		       << " will not be constructed \n";
	          throw NoEnergiesDefined(diag.str());

              }
           }
        }
        ++s;        
    }     

    int Nchannels (0);
    size_t numModVects = m_modelPositionDirectory[groupIndex-1].size();
    if (groupSpectra.empty()) 
    {
        m_groups.push_back(new PlotGroup(groupIndex,0, numModVects) );
    }
    else
    {

        Nchannels = (*groupSpectra.begin())->channels();      
        // create a PlotGroup & allocate its arrays.
        m_groups.push_back(new PlotGroup(groupIndex,Nchannels,numModVects) ); 
    }

    PlotGroup* gr (m_groups.back());
    gr->single = (Nchannels == 1);

    std::map<int, RebinInfo>::iterator rbEntry = m_groupsRebinInfo.find(groupIndex); 
    if (rbEntry == m_groupsRebinInfo.end())
    {
	// Use all-group default values, key = -1, always in map.
	rbEntry = m_groupsRebinInfo.find(-1);
    }
    gr->critSignificance = rbEntry->second.sigma;
    gr->critNumChans = rbEntry->second.maxBins;
    gr->errorType    = rbEntry->second.mode; 

    std::list<SpectralData*>::const_iterator gs = groupSpectra.begin();
    std::list<SpectralData*>::const_iterator gsEnd = groupSpectra.end();

    int maxChan (0);
    int minChan (999999);


    while ( gs != gsEnd )
    {
        int j (0);
        const BoolArray& nc = (*gs)->noticedChannels();        

        BoolArray::const_iterator it = nc.begin();
        BoolArray::const_iterator Nend = nc.end();
        while ( !nc[j] && it != Nend ) ++it, ++j;     
        minChan = std::min(minChan, j);

        it = nc.end();
        j  = nc.size();
        do
        {
	    --it, --j;                                
        }       while (!nc[j] && it != nc.begin()); 

        maxChan = std::max(maxChan,j);
        ++gs;       
    }

    // chanXSpacing is used to separate plot groups when they are
    // viewed in channel space.

    gr->channelPlotSpacing = (maxChan - minChan);
    if ( gr->channelPlotSpacing <= 0 )
    {
        // group is totally ignored.
        gr->channelPlotSpacing = 10;
        gr->lastInputChannel = 1;
    }
    else if ( gr->single )
    {
        // single point (data contains observation through a filter,
        // not a spectrum).
        gr->channelPlotSpacing = 20;
        gr->lastInputChannel = 1;
    }
    else
    {
        // square up to nearest multiple of 10.
        gr->channelPlotSpacing = 10*( (gr->channelPlotSpacing+2)/10 + 1 );       
        gr->lastInputChannel = maxChan + 1;
    }

    const std::map<size_t,size_t>& sourceAssignments = 
                m_modelPositionDirectory[gr->plotGroup - 1];
    std::map<size_t,size_t>::const_iterator itAssignment =
                sourceAssignments.begin();
    std::map<size_t,size_t>::const_iterator itAssignmentEnd =
                sourceAssignments.end();
    while (itAssignment != itAssignmentEnd)
    {
       size_t sourceNum = itAssignment->first;
       size_t sourcePos = itAssignment->second;
       // Any info in sourceAssigments has already been verified to
       // correspond to an active model object.  No need to check things
       // again in here.
       const string modName = models->lookupModelForSource(sourceNum);
       const Model* mod = models->lookup(modName);
       const size_t compOffset = gr->model.size() > 1 ? 1 : 0;
       gr->model[sourcePos].lineStyle = m_modelStyle;

       if (m_showAddComponents || unfold)
       {
          size_t numSumComps = mod->sources().size(); 
          if (numSumComps > 1)
          {
             PlotVector* x = &(gr->xAxis);
             std::list<PlotVector>& currentComps =
                   gr->sources[sourcePos-compOffset];
             for (size_t i=0; i<numSumComps; ++i)
             {
                currentComps.push_back(PlotVector(gr->n));
                currentComps.back().xData = x;
             }
	     // If unfold, need to allocate arrays for the unfolded
	     // versions of the source components.  
	     if ( unfold )
	     {
	         std::list<PlotVector>& currentCompFluxes 
		     = gr->sourceFluxes[sourcePos-compOffset];
	         for (size_t j = 0; j < numSumComps; ++j)
	         {
		     currentCompFluxes.push_back(PlotVector(gr->n)); 
		     currentCompFluxes.back().xData = x;
	         }                                              
	     }
          }
       } // end if show add components      
       ++itAssignment;
    }
    // Sum of multiple models will be placed in the 0th element.
    if (sourceAssignments.size() > 1)
       gr->model[0].lineStyle = m_modelStyle;

    return gr;
}

void Plot::integrateModelFluxes (std::list<SpectralData*>& spectra, Real eMin, Real eMax, PlotGroup* group, int point)
{
    // encapsulate the following in a private member function
    // integrate flux
    const std::map<size_t,size_t>& modAssignments =
          m_modelPositionDirectory[group->plotGroup-1];
    const size_t compOffset = group->model.size() > 1 ? 1 : 0;
    std::list<SpectralData*>::const_iterator gs    (spectra.begin());
    std::list<SpectralData*>::const_iterator gsEnd (spectra.end());     
  // if there is more than one source per spectrum, the model fluxes are stored
  // in the flux std:vector<Real*> container enumerated from 1.
  // if there isn't then that vector only has one element, enumerated as zero.
  // for the latter case, the source number is 1 and the binned flux is stored
  // in location zero. for the former case, the binned fluxes are stores in
  // location == source number and the sum will ultimately be stored in location
  // zero.

    while (gs != gsEnd)
    {
       size_t spectrumNumber = (*gs)->spectrumNumber();
       std::vector<const Model*>::const_iterator itMod = 
            m_modelsForSpectra[spectrumNumber-1].begin();
       std::vector<const Model*>::const_iterator itModEnd = 
            m_modelsForSpectra[spectrumNumber-1].end();
       while (itMod != itModEnd)
       {
          const Model* m = *itMod;
          const size_t sourceIndex = 
                    modAssignments.find(m->sourceNumber())->second;
          const ArrayContainer& flux = m->modelFlux();
          if ( flux.find(spectrumNumber) != flux.end())
          {
             Real energyFlux(0);
             Real ergFlux(0);
             m->integrateFlux(spectrumNumber,eMin,eMax,energyFlux,ergFlux);
             Real width(1.);
             if ( std::abs((eMax - eMin)/eMax) > 1.E-15)
             {
                switch (m_xOption)
                {
                   case CHANNELS:
                   case ENERGY:
                   default:
                           width /= (eMax - eMin);
                   break;          
                   case WAVELENGTH:
                      if (m_waveUnitsHz == 1)
                         width /= Numerics::KEVTOHZ*(eMax - eMin);
                      else
                         width /= Numerics::KEVTOA*(1./eMin - 1./eMax);
                   break;
                }
                group->auxData[sourceIndex].data[point] 
                                = energyFlux*width;
             }
             else
             {
                group->auxData[sourceIndex].data[point] 
                                     =  s_NODATA;
             }

             // compute the integrals also for the source components, if there are any.
             if (!group->sourceFluxes[sourceIndex-compOffset].empty() )
             {
                 std::list< std::pair<Real,Real> > fluxes 
                     = m->integrateSourceFlux(spectrumNumber,eMin,eMax); 
                 std::list<std::pair<Real, Real> >::const_iterator fl = fluxes.begin();
                 std::list<std::pair<Real, Real> >::const_iterator flEnd = fluxes.end();
                 PlotVectorList::iterator sfl = group->sourceFluxes[sourceIndex-compOffset].begin();
                 if ( std::abs(eMax - eMin) > 1.E-5)
                 {                            
                    while ( fl != flEnd )
                    {
                       sfl->data[point] = (fl->first)*width;
                       ++fl, ++sfl;       
                    }
                 }
                 else
                 {
                    while ( fl != flEnd )
                    {
                       sfl->data[point] = s_NODATA;
                       ++fl, ++sfl;       
                    }                                        
                 }
             }

          }
          ++itMod;
       }                                                            
       ++gs;
    }
}

void Plot::setXaxis (PlotGroup* group, int firstChannel, int lastChannel, Real eMin, Real eMax, int point)
{

    switch (m_xOption)
    {
          default:
          case CHANNELS:
                if ( group->single )
                {
                        group->xAxis.errors[0][point] = 0.5*group->channelPlotSpacing;      
                        group->xAxis.data[point] = group->xAxis.errors[0][point]  - 0.5;
                        group->xAxis.data[point] = group->xAxis.data[point];            
                }
                else
                {
                        group->xAxis.data[point] =  0.5*(lastChannel + firstChannel); 
                        group->xAxis.data[point] = group->xAxis.data[point];            
                        group->xAxis.errors[0][point] =  0.5*(lastChannel - firstChannel);                       
                }
                break;
          case ENERGY:
                group->xAxis.data[point] = 0.5*(eMax + eMin);
                group->xAxis.errors[0][point] = 0.5*(eMax - eMin);
                break;
          case  WAVELENGTH:
                if ( eMax <= 0 || eMin <= 0)
                {
                        group->xAxis.data[point] = s_NODATA;
                        group->xAxis.errors[0][point] = s_NODATA;       
                }
                else
                {
                        group->xAxis.data[point]      = 0.5*Numerics::KEVTOA*(1./eMax + 1./eMin);
                        group->xAxis.errors[0][point] =  0.5*Numerics::KEVTOA*(1./eMin - 1./eMax);                                     
                }
                break;  
    }  
}

void Plot::setError (Plot::PlotErrMode errorType, Real& error, Real& backError, Real areaSum, Real backAreaSum, Real count, Real backgroundCount, Real bscaleRatio)
{
  switch ( errorType )
  {
        case ROOTN:
                error = sqrt(count + backgroundCount);
                backError = sqrt(backgroundCount);
                if (!m_counts) 
                {
                   error /= areaSum;
                   backError /= areaSum;
                }
                break;
        case GEHRELS1:
                error = (1 + sqrt(0.75+count))*(1 + sqrt(0.75+count));
                if ( backgroundCount != 0) 
                {
                   backError = (1 + sqrt(0.75+backgroundCount))*(1 + sqrt(0.75+backgroundCount));      
                   error += backError;
                }
                else 
                   backError = 0.0;
                if (!m_counts)
                {
                   error /= areaSum;
                   backError /= areaSum; 
                }               
                break;
        case GEHRELS2:
                error = count - 0.25;
                if ( backgroundCount != 0)
                {
                   backError = backgroundCount - 0.25;
                   error += backError;
                }
                else
                   backError = 0.0;
                error = sqrt(error);
                backError = sqrt(backError);
                if (!m_counts)
                {
                   error /= areaSum; 
                   backError /= areaSum;
                }   
                break;
        case GEHRELSM:
                error = ((1 + sqrt(count+0.75) + sqrt(count - 0.25))*0.5)*
                          ((1 + sqrt(count+0.75) + sqrt(count - 0.25))*0.5);
                if ( backgroundCount != 0)
                {
                   backError = ((1 + sqrt(backgroundCount + 0.75) 
                                   + sqrt(backgroundCount - 0.25))*0.5)*
                            ((1 + sqrt(backgroundCount + 0.75) 
                                   + sqrt(backgroundCount - 0.25))*0.5);
                   error += backError;                                               
                }
                else
                   backError = 0.0;  
                error = sqrt(error);
                backError = sqrt(backError);
                if (!m_counts)
                {
                   error /= areaSum; 
                   backError /= areaSum;
                }   
                break;          
        case STD:
        default:
                error = sqrt(error);
                backError = sqrt(backError);
                break;

  }

  // If chisq is the current statistic, plot chisq and delchi can't
  // have zero errors.
  if (dynamic_cast<ChiSquare*>(XSContainer::fit->statMethod()))
  {
     if (std::fabs(error) <= SMALL)
     {
         error = ChiSquare::calcMinVariance(areaSum, backAreaSum, bscaleRatio);
         error = sqrt(error);
         if (m_counts)
            error *= areaSum;
     }
  }
}

void Plot::setSelectedLabel (size_t nLabel, const string& labelKey, const string& option)
{
  if (nLabel >= m_selectedLabels.size())
  {
     // nLabel is not user input.  This situation should
     // never exist.  If it does, something is seriously wrong.
     throw RedAlert("Error loading plot labels");
  }

  if (labelKey.size())
  {
     string label(labelNames(labelKey));
     if ( option.length())
     {
           label += " ";
           label += option;
     }
     // This throws a YellowAlert if labelKey is not in map.
     m_selectedLabels[nLabel] = label ;
  }
  else
  {
     // Allow option of deleting a label if key is empty.
     m_selectedLabels[nLabel] = "";
  }
}

void Plot::fixYMin (Real value)
{
  m_ranges["y"].first = value;
}

void Plot::autoSelectedLabels (const string& key)
{

  // Automatic generation of x, y, based on the command dependent
  // key input.  To be used only when the x axis label is dependent 
  // on the XOption setting. An asterisk at the end of the key indicates 
  // that the y label is also dependent on the current XOption setting. 

  // xAxis label sets the x-axis label and returns the character
  // representing the axis type (channels, energy or wave).
  char xchar(xAxisLabel());      

  // Added code: Allow normalized by area option: if key is
  // is appended by "(whitespace)a".
  string subKey = key;
  string inverseArea;
  string::size_type aLoc = key.find(" a");
  if (aLoc != string::npos)
  {
     // remove " a" from the key
     subKey = key.substr(0,aLoc);
     inverseArea = labelNames("normArea"); 
  }

  string fullKey ("y" + subKey);
  size_t i = fullKey.length();
  if (fullKey[i-1] == '*')
  {
     fullKey[i-1] = xchar;
  }
  // Another addition: If xoption is wave, decide between Hz units
  // (default) and A units.
  if (fullKey[i-1] == 'w')
  {
     if (m_waveUnitsHz == 0)
     {
        fullKey += 'A';
     }
  }

  setSelectedLabel(1, fullKey, inverseArea);
}

void Plot::initializeSpectra (XSContainer::DataContainer* dataContainer)
{

  // all non-throwing operations - except for the RedAlert which,
  // of course, is a terminate call.      

  m_spectra.clear();

  size_t N (dataContainer->numberOfSpectra());
  for ( size_t j = 1; j <= N; ++j)
  {
        SpectralData* sp (datasets->lookup(j));       
        if (sp)
        {
                m_spectra.insert(SpecGroup::value_type(sp->plotGroup(),sp));
        }
        else
        {
                throw RedAlert(" data container is corrupted: missing spectra ");
        }   
  }

}

void Plot::resetPlotGroups ()
{

  PlotGroupContainer::iterator doomed = m_groups.begin();
  PlotGroupContainer::iterator gEnd = m_groups.end();

  while ( doomed != gEnd )
  {
        delete *doomed;
        ++doomed;    
  }

  m_groups.clear();
  m_plotBeginIndex = 0;
}

void Plot::resetAttributes ()
{
  // This function is intended to reset only those attributes of
  // Plot that the Command objects shouldn't have to concern themselves
  // with initializing everytime they are called.  Things the command
  // objects are expected to set, such as the default range flags, should
  // not be set here.

  m_xLine = std::pair<bool,Real>(false,.0);
  m_numGraphs = 1;
  m_currentGraphNumber = 0;
}

void Plot::bundlePlotVectors (bool addComponents)
{
   PlotGroupContainer::iterator iPg = m_groups.begin();
   PlotGroupContainer::iterator iPgEnd = m_groups.end();
   size_t iGroup=0;
   Real yLow = 1.0e32;
   Real yHigh = -1.0e32;
   while (iPg != iPgEnd)
   {
      if (iGroup >= m_plotBeginIndex)
      {
         PlotGroup* gr = (*iPg);
         std::vector<PlotVector*>& groupAttributes = gr->bundledPlotVectors;
         groupAttributes.clear();
         groupAttributes.push_back(&gr->xAxis);
         int modStart = 1;
         if (gr->yData.lineStyle || gr->yData.symbolStyle)
         {
            groupAttributes.push_back(&gr->yData);
            ++modStart;
         }
         if (gr->background.lineStyle || gr->background.symbolStyle)
         {
            groupAttributes.push_back(&gr->background);
            ++modStart;
         }
         gr->plotVectorBoundaries[0] = modStart;
         int compStart = modStart;
         for (size_t j=0; j < gr->model.size(); ++j)
         {
            if (gr->model[j].lineStyle || gr->model[j].symbolStyle)
            {
               groupAttributes.push_back(&gr->model[j]);
               ++compStart;
            }
         }
         gr->plotVectorBoundaries[1] = compStart;
         int auxStart = compStart;      

         if (m_showAddComponents || addComponents ) 
         {
            for (size_t j=0; j < gr->sources.size(); ++j)
            {
                 PlotVectorList::iterator sl (gr->sources[j].begin());
                 PlotVectorList::iterator slEnd (gr->sources[j].end());
                 while ( sl != slEnd )
                 {
                    if (sl->lineStyle || sl->symbolStyle)
                    {
                       groupAttributes.push_back( &(*sl) );
                       ++auxStart; 
                    }                 
                    ++sl;
                 }
             }
         }  

         gr->plotVectorBoundaries[2] = auxStart;    
         for (size_t j=0; j<gr->auxData.size(); ++j)
         {
            if (gr->auxData[j].lineStyle || gr->auxData[j].symbolStyle)
               groupAttributes.push_back(&gr->auxData[j]);
         }

         // The setting of the y range for only the top graph is a
         // non-interface-changing patch hack to implement setplot id.
         // (See PlotCommand and Plt for the rest of it.)  This is 
         // to determine the y-limits for the commands which rely
         // on getNextPlotBin.  For the other commands this is
         // currently redundant, since they already set a y range
         // in their command handlers.  Still, redoing it here
         // shouldn't hurt anything.         
         if (m_currentGraphNumber == 0 && m_plotIDs &&
                m_xOption != Plot::CHANNELS)
         {
            for (size_t j=1; j<groupAttributes.size(); ++j)
            {
               typedef std::vector<Real>::const_iterator YIter;
               YIter itY = groupAttributes[j]->data.begin();
               YIter itYEnd = groupAttributes[j]->data.end();
               YIter itMax = std::max_element(itY, itYEnd);
               YIter itMin = std::min_element(itY, itYEnd);
               if (*itMax > yHigh)   yHigh = *itMax;
               if (*itMin < yLow)   yLow = *itMin;
            }
            m_ranges["y"].first = yLow;
            m_ranges["y"].second = yHigh;
         }

      }
      ++iGroup;
      ++iPg;
   }
}

char Plot::xAxisLabel ()
{
  char xchar;  
  switch (m_xOption)
  {
     default:
     case CHANNELS:
	xchar = 'c';
	break;
     case ENERGY:
	xchar = 'e';
	break;
     case WAVELENGTH:
	xchar = 'w';
	break;
  }
  string fullKey = "x";
  fullKey += xchar;
  setSelectedLabel(0, fullKey);
  return xchar;
}

Plot::ModelStatus Plot::determineModelStatus ()
{
  const ModelMap& mods = models->modelSet();
  ModelMapConstIter mm (mods.begin());
  ModelMapConstIter mmEnd (mods.end());

  m_modelStatus = INACTIVE;
  while (mm != mmEnd)
  {
     if (mm->second->isActive())
     {
        m_modelStatus = FOLDED;
	break;
     }
     ++mm;
  }
  return m_modelStatus;
}

void Plot::makeTheGraph ()
{
  try
  {
     m_graph->perform(this);
  }
  catch (YellowAlert &)
  {
     m_grid = 0;
     m_savedPlotArray.clear();
     resetAttributes();
     resetPlotGroups();
     m_graph->flushHardcopy();
     throw;
  }
  m_grid = 0; // non-owning pointer
  saveArray();
  resetAttributes();
  resetPlotGroups();
  m_graph->flushHardcopy();
}

void Plot::prepareTheGraph ()
{
   m_graph->prepare(this);
}

PlotGroup* Plot::initializePlotGroup (SpectralData* spectrum)
{
  // Version of initializePlotGroup for handling efficiency, sensitivity, and
  // insensitivity plot commands.
  PlotGroup* gr(0); 
  int spectrumNumber (spectrum->spectrumNumber());        
  bool increasing ( m_xOption == CHANNELS || m_xOption == ENERGY  );  
  const std::vector<Response*>& dets = spectrum->detector();
  bool anyDets = false;
  for (size_t iDet=0; iDet<dets.size(); ++iDet)
  {
     const Response* resp = dets[iDet];
     if (resp)
     {
        anyDets = true;
        const RealArray& energy = resp->energies();
        size_t NE (energy.size() - 1);
        m_groups.push_back(new PlotGroup(m_groups.size()+1,NE,1));
        gr = m_groups.back();

        for (size_t j = 0; j < NE; ++j)
        {
              // dummies are only used for channel mode, which has to be
              // switched off.
              int dummy(0);
              setXaxis(gr,dummy,dummy,energy[j],energy[j+1],j);      
        }      

        if ( increasing )
        {
                m_ranges["x"].first = std::min(gr->xAxis.data[0] - gr->xAxis.errors[0][0],
                                      m_ranges["x"].first);
                m_ranges["x"].second = std::max(gr->xAxis.data[NE-1] + gr->xAxis.errors[0][NE-1],
                                      m_ranges["x"].second);
        }
        else
        {
                m_ranges["x"].first = std::max(gr->xAxis.data[NE-1] + gr->xAxis.errors[0][NE-1],
                                      m_ranges["x"].first);
                m_ranges["x"].second = std::min(gr->xAxis.data[0] - gr->xAxis.errors[0][0],
                                      m_ranges["x"].second);
        }

     } // end if resp
  } // end det loop

  if (!anyDets)
  {
        // one could put in a dummy response here, but for the uses
        // this function is likely to be put to the plot won't be very
        // interesting (e.g. detector efficiency).
        // The Exception is 
        std::ostringstream ss;
        ss << spectrumNumber << " response not loaded \n";
        throw NoEnergiesDefined(ss.str());
  }

  return gr;
}

PlotGroup* Plot::initializePlotGroup (Model* model, bool addComponents)
{
  PlotGroup* gr (0);
  const ArrayContainer& energies = model->energy();

  bool increasing ( m_xOption == CHANNELS || m_xOption == ENERGY  );

  ArrayContainer::const_iterator en = energies.begin();
  ArrayContainer::const_iterator enEnd = energies.end();

  while ( en != enEnd )
  {
     int spectrumNumber (en->first);
     const RealArray& energy = en->second;
     int NE (energy.size() - 1);    
     m_groups.push_back(new PlotGroup(m_groups.size()+1,NE, 1));
     gr = m_groups.back();

     // set the x axis array.

     for (int j = 0; j < NE; ++j)
     {
        // dummies are only used for channel mode, which has to be
        // switched off.
        int dummy(0);
        setXaxis(gr,dummy,dummy,energy[j],energy[j+1],j);      
     }      

     if ( increasing )
     {
        m_ranges["x"].first  = std::min(gr->xAxis.data[0] 
                 - gr->xAxis.errors[0][0], m_ranges["x"].first);
        m_ranges["x"].second = std::max(gr->xAxis.data[NE-1] 
                 + gr->xAxis.errors[0][NE-1], m_ranges["x"].second);
     }
     else
     {
        m_ranges["x"].second = std::min(gr->xAxis.data[0] 
                 - gr->xAxis.errors[0][0], m_ranges["x"].second);
        m_ranges["x"].first  = std::max(gr->xAxis.data[NE-1] 
                 + gr->xAxis.errors[0][NE-1], m_ranges["x"].first);
     }

     const RealArray& flux = model->modelFlux(spectrumNumber);
     const RealArray& fluxErr = model->modelFluxError(spectrumNumber);

     gr->model[0].data.resize(NE,0);
     gr->model[0].lineStyle = PlotStyle::SOLID;
     for ( int j = 0; j < NE; ++j)
     {
        gr->model[0].data[j] = flux[j];
     }
     if ( fluxErr.size())
     {
        gr->model[0].errors.push_back(std::vector<Real>(NE));
        for ( int j = 0; j < NE; ++j)
        {
           gr->model[0].errors[0][j] = flux[j];
        }     
     }
     else 
     {
        gr->model[0].errors.push_back(std::vector<Real>());             
     }

     if ( m_showAddComponents || addComponents )
     {
        // if sources are present and required, put them in the right 
        // place too.

        const SourceList& sources = model->sources();
        size_t numberOfSourcesInModel(sources.size()); 
        if (numberOfSourcesInModel > 1)
        {
           const size_t& n = gr->n;
           PlotVector* x = &(gr->xAxis);
	   PlotVectorList& currentSources = gr->sources[0];

           SourceList::const_iterator sl = sources.begin();    

           for (size_t j = 0; j < numberOfSourcesInModel; ++j)
           {
               currentSources.push_back(PlotVector(n));
               PlotVector& cs = currentSources.back();
               const SumComponent* c = *sl;
               for (size_t k = 0; k < n; ++k)
               {
                  cs.data[k]  = c->photonArray(spectrumNumber)[k];       
               }
               cs.xData = x;
               cs.lineStyle = PlotStyle::DASHDOT;
               ++sl;
           }
        }
     }
     ++en;  
  }

  return gr;
}

void Plot::makePlotArraysFromModels (const string& name, bool addComponents)
{
  // unbinned plot arrays makes one PlotGroup object per spectrum
  // and ignores the plot grouping mechanism.
  if ( m_xOption == CHANNELS || m_xOption == ENERGY )
  {
       m_ranges["x"].first  =  1.E+100;    
       m_ranges["x"].second = -1;    
  }
  else
  {
       m_ranges["x"].first  = -1;    
       m_ranges["x"].second = 1.E+100;             
  }  

  resetPlotGroups();

  string mapName(name);
  if ( !mapName.size() ) mapName = Model::DEFAULT();

  if (models->modelSet().empty())
  {
     throw NoModelDefined("");
  }
  std::pair<ModelMapConstIter, ModelMapConstIter> im (models->modelSet().equal_range(name));

  if (im.first == im.second)
  {
     string msg;
     if (name == Model::DEFAULT())
     {
        msg = "No model has been entered with the default name";
     }
     else
     {
        msg = "No model has been entered with name = " + name;
     }
     throw ModelNotFound(msg);
  }  

  if (datasets->numSourcesForSpectra() && im.first->second->sourceNumber() > datasets->numSourcesForSpectra())
  {
     string msg("Cannot plot model of higher source number than data allows.");
     throw NotImplemented(msg);
  }

  m_ranges["y"].first  = LARGE;
  m_ranges["y"].second = 0;

  ModelMapConstIter m (im.first);

  while ( m  != im.second )
  {

        initializePlotGroup(m->second, addComponents);
        ++m;

  }
}

void Plot::fixXMin (Real value)
{
  m_ranges["x"].first = value;
}

void Plot::fixXMax (Real value)
{
  m_ranges["x"].second = value;
}

void Plot::fixYMax (Real value)
{
  m_ranges["y"].second = value;
}

void Plot::setRangesFromGrid ()
{

  // will actually throw if grid is not yet initialized before this routine
  // is called, but to get a total safety guarantee, add the 'if'.
  if ( m_grid )
  {
        const Grid::SpecContainer& paramArray = m_grid->getParameter();
        // ... then Ndim must be >= 1 .
        Grid::SpecContainer::const_iterator p = paramArray.begin();
        Grid::ParameterSpec* parameter1  (*p);
        // unroll the first three... (?) don't really know what to do with a 4-D grid
        // and can't plot a 3-D now anyway, but...
        // NOTE: Ranges are NOT taken from the lowRange and highRange
        // members of ParameterSpec.  Those correspond to the user-
        // entered ranges, which may then be shifted to bin centers
        // (ie. as with margin).  The parameterValues array takes
        // this into account.
        m_ranges["x"].first  = parameter1->parameterValues.front();
        m_ranges["x"].second = parameter1->parameterValues.back();
        ++p;
        if ( p == paramArray.end() )
        {
           const RealArray& values = m_grid->getGrid();
           m_ranges["y"].first  = values.min();
           m_ranges["y"].second = values.max();                
        }
        else
        {
           Grid::ParameterSpec* parameter2 (*p);
           m_ranges["y"].first  = parameter2->parameterValues.front();
           m_ranges["y"].second = parameter2->parameterValues.back();   
           ++p;
           if ( p != paramArray.end() )
           {
              Grid::ParameterSpec* parameter3 (*p);
              m_ranges["z"].first  = parameter3->parameterValues.front();
              m_ranges["z"].second = parameter3->parameterValues.back();   
              ++p;
              int h(1);
              while ( p != paramArray.end() )
              {
                 std::ostringstream dimName;
                 dimName << "t" << h;
                 m_ranges[dimName.str()].first = (*p)->parameterValues.front();       
                 m_ranges[dimName.str()].second = (*p)->parameterValues.back();       
                 ++p, ++h;
              }        
           }
        }
  }      
}

void Plot::ungroupAll ()
{
  int sz = datasets->numberOfSpectra();
  for (int i=0; i<sz; ++i)
  {
     datasets->setPlotGroupNums(i, i+1);
  }
  // This call is needed to update datasets' m_numberOfPlotGroups.
  datasets->renumberPlotGroups(1, 0);
  initializeSpectra(datasets);
}

void Plot::initEffectiveAreas (std::list<SpectralData*>& spectraInGroup)
{
   std::list<SpectralData*>::iterator it = spectraInGroup.begin();
   std::list<SpectralData*>::iterator itEnd = spectraInGroup.end();
   while (it != itEnd)
   {
      (*it)->calcEffectiveAreas();
      ++it;
   }
}

void Plot::clearEffectiveAreas (std::list<SpectralData*>& spectraInGroup)
{
   std::list<SpectralData*>::iterator it = spectraInGroup.begin();
   std::list<SpectralData*>::iterator itEnd = spectraInGroup.end();
   while (it != itEnd)
   {
      (*it)->clearEffectiveAreas();
      ++it;
   }
}

void Plot::saveArray ()
{
   // Do not throw any YellowAlerts from here.  If saveArrayInfo
   // does not correspond to an existing plot group or array, 
   // just issue warning and leave m_savedPlotArray empty.
   if (m_savedPlotArray.size())  m_savedPlotArray.resize(0);
   if (m_saveArrayInfo.first != NO_SAVE)
   {
      int groupNum = m_saveArrayInfo.second;
      if (groupNum < 1 || static_cast<size_t>(groupNum) > m_groups.size())
      {
         tcout << "***Warning: Saved plot array group number specifier = "
           << groupNum << ",\n   does not correspond to an existing plot group.  "
           << "No plot array will be saved." << std::endl; 
      }
      else
      {
         const PlotGroup* pg = m_groups[groupNum - 1];
         if (pg->plotGroup != groupNum)
         {
            throw RedAlert("In Plot::saveArray: plot group nums not sequentially numbered.");
         }
         switch (m_saveArrayInfo.first)
         {
            case X_AR:
               m_savedPlotArray = pg->xAxis.data;
               break;
            case XERR_AR:
               m_savedPlotArray = pg->xAxis.errors[0];
               break;
            case Y_AR:
               m_savedPlotArray = pg->yData.data;
               break;
            case YERR_AR:
               m_savedPlotArray = pg->yData.errors[0];
               break;
            case MODEL_AR:
               // I doubt we can get here with an empty vector of 
               // model PlotVectors.  Still to be safe...
               if (pg->model.size())
               {
                  m_savedPlotArray = pg->model[0].data;
               }
               else
               {
                  tcout << "***Warning: No model plot array exists." << std::endl;
               }
               break;
            default:
               break;
         }
      }
   }
}

void Plot::itemizePlotGroupModels ()
{
   // Only count active/on models.
   m_modelPositionDirectory.clear();
   m_modelsForSpectra.clear();
   // This is the highest possible number of models used.
   const size_t allSources = datasets->numSourcesForSpectra();
   // To speed things along so we don't always have to check every detector
   // slot for every spectrum, first see which slots actually have active 
   // models associated with them, and store the number of sum comps for 
   // a given model so its lookup only needs to be performed once.
   std::vector<size_t> sourcesWithActiveMods;
   std::vector<size_t> numCompsForMods(allSources,0);
   for (size_t i=0; i<allSources; ++i)
   {
      string modName(models->lookupModelForSource(i+1));
      // Model with modName must be active/on or active/off.
      if (modName.length())
      {
         const Model* mod = models->lookup(modName);
         if (mod && mod->isActive())
         {
            // Let's keep source numbers 1-based throughout.
            sourcesWithActiveMods.push_back(i+1);
            const size_t nComps = mod->sources().size();
            numCompsForMods[i] = nComps;
         }
      }
   }
   size_t prevGroup = 1;
   m_modelsForSpectra.resize(m_spectra.size());
   std::set<size_t> sourcesForGroup;
   SpecGroup::const_iterator itSpec = m_spectra.begin();
   SpecGroup::const_iterator itSpecEnd = m_spectra.end();
   while (itSpec != itSpecEnd)
   {
      // Don't assume spectra are in order by spectrum number,
      // only by plot group number.
      const std::vector<const Model*> 
                modsForSpec(models->getModsForSpec(itSpec->second));
      std::vector<const Model*>& spectraMods = 
                m_modelsForSpectra[itSpec->second->spectrumNumber()-1];
      // modsForSpec has by definition size = allSources, and
      // may have many (or all) elements = 0.  We don't need to
      // check each one, thanks to sourcesWithActiveMods subset.
      std::vector<size_t>::const_iterator itActSource =
                sourcesWithActiveMods.begin();
      std::vector<size_t>::const_iterator itActSourceEnd =
                sourcesWithActiveMods.end();
      while (itActSource != itActSourceEnd)
      {
         const Model* mod = modsForSpec[*itActSource-1];
         if (mod)
         {
            spectraMods.push_back(mod);
            // Don't want to repeat source nums, that's
            // why this is a set and not a multiset.
            sourcesForGroup.insert(mod->sourceNumber());
         }
         ++itActSource;
      }
      ++itSpec;
      if (itSpec == itSpecEnd || itSpec->first != prevGroup)
      {
         // End of current plot group.
         if (itSpec != itSpecEnd)
            prevGroup = itSpec->first;
         // If more than one model, reserve 0 element for sum.
         size_t offset = sourcesForGroup.size() > 1 ? 1 : 0;
         std::set<size_t>::const_iterator itNum = sourcesForGroup.begin();
         std::map<size_t,size_t> assignments;
         size_t count = 0;
         while (itNum != sourcesForGroup.end())
         {
            size_t sourceForGroup = *itNum;
            assignments[sourceForGroup] = count + offset;
            ++itNum;
            ++count;
         }
         m_modelPositionDirectory.push_back(assignments);
         sourcesForGroup.clear();
      }
   }
}

// Additional Declarations