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

// DataContainer
#include <XSModel/GlobalContainer/DataContainer.h>
// ModelContainer
#include <XSModel/GlobalContainer/ModelContainer.h>
// PlotInsensitivity
#include <XSPlot/Commands/PlotInsensitivity.h>
#include <XSPlot/Plot/Plot.h>
#include <XSModel/Data/SpectralData.h>
#include "XSContainer.h"
#include "XSstreams.h"
#include "XSsymbol.h"


// Class PlotInsensitivity 

PlotInsensitivity::PlotInsensitivity (Plot* plot)
  : PlotCommand(plot)
{
    modelRequired(Plot::ACTIVE);
}


void PlotInsensitivity::manipulate ()
{
  using namespace XSContainer;
  static const Real PHOT_TO_MJY (662.22);
  static const Real BIG (1.E10);

  Plot* pPlot = plot();
  pPlot->ranges("y").first  = LARGE;
  pPlot->ranges("y").second  = -LARGE;

  PlotGroupContainer::iterator itPg ( pPlot->groups().begin());
  PlotGroupContainer::iterator itPgEnd ( pPlot->groups().end());
  bool anyModels = false;
  const size_t NS (datasets->numberOfSpectra());
  for (size_t iSpec=1; iSpec<=NS; ++iSpec)
  {
     // A separate plot group must exist for every response of every spectrum,
     // though not all of these necessarily have an active model associated
     // with them. If they don't, they will be turned off.
     SpectralData* spectrum = datasets->lookup(iSpec);
     std::vector<const Model*> specMods = models->getModsForSpec(spectrum);
     const std::vector<Response*>& detectors = spectrum->detector();
     const size_t N = detectors.size();
     std::vector<RealArray> sensitivity(spectrum->sensitivity());
     for (size_t iDet=0; iDet<N; ++iDet)
     {
        size_t l=0;
        if (detectors[iDet])
        {
           PlotGroup* gr(*itPg);
           RealArray& thisSensitivity = sensitivity[l];
           const RealArray& energy = 
                const_cast<const Response*>(detectors[iDet])->energies();
           if (specMods[iDet])
           {
              anyModels = true;
              for (size_t k=0; k<gr->n; ++k)
              {
                 Real ins(PHOT_TO_MJY*sqrt(energy[k]*energy[k+1]/thisSensitivity[k]));
                 if (ins < BIG)
                 {
                    gr->model[0].data[k] = ins;
                    pPlot->ranges("y").first 
                               = std::min(gr->model[0].data[k],
                                    plot()->ranges("y").first);
                    pPlot->ranges("y").second 
                               = std::max(gr->model[0].data[k],
                                    plot()->ranges("y").second);
                 }
                 else
                 {
                    gr->model[0].data[k] = BIG;
                 }
              }
              gr->model[0].lineStyle = PlotStyle::DASHDOT;       
           }
           else
           {
              gr->model[0].lineStyle = PlotStyle::NONE;       
           }
           gr->model[0].symbolStyle = PlotStyle::BLANK;       
           gr->yData.lineStyle = PlotStyle::NONE;
           gr->yData.errors[0].resize(0);


           ++itPg, ++l;
        }
     } 
  }
  if (!anyModels)
  {
     throw YellowAlert("\nCannot plot sensitivity, no models found.\n");
  }

  // Anything that was set to BIG should now be reset to the max y val.
  Real maxYval = pPlot->ranges("y").second;
  itPg = pPlot->groups().begin();
  while (itPg != itPgEnd)
  {
     std::vector<Real>::iterator itVal = (*itPg)->model[0].data.begin();
     std::vector<Real>::iterator itValEnd = (*itPg)->model[0].data.end();
     while (itVal != itValEnd)
     {
        Real& val = *itVal;
        if (val > maxYval)
           val = maxYval;
        ++itVal;
     }
     ++itPg;
  }
}

void PlotInsensitivity::primaryGraph (const string& value)
{
  Plot* pPlot = plot();
  Plot::XaxisMode saveXoption (pPlot->xOption());
  pPlot->setXOption("energy");
  if (saveXoption == Plot::CHANNELS)
  {
     pPlot->plotXLog(pPlot->xLog());
  }
// set titles
  pPlot->setSelectedLabel(0,"xe");
  pPlot->setSelectedLabel(1,"yy");
  pPlot->setSelectedLabel(2,"ty");

  try
  {
   // manufacture necessary data arrays into pPlot's PlotGroup container.
   // For this particular graph need X axis set by setXaxis using the 
   // Response's energy array (Ehi, Elo, not EboundsMax, EboundsMin)
   // and its efficiency array.

     pPlot->makePlotArraysFromSpectra();

   // the manipulate call manipulates the structures produced by makeUnbinnedPlotArrays
   // to get the dedicated arrays required by this particular plot. In this case, it
   // grabs the efficiency arrays from the Response.

     manipulate();
     fillModelAttributes(1,1);
     pPlot->fixYMax(pPlot->ranges("y").second*1.01);
     pPlot->rangeDefaultFlags(3);

   // graph it.

     pPlot->bundlePlotVectors();
     pPlot->prepareTheGraph();
     pPlot->makeTheGraph(); 
  }
  catch (YellowAlert &)
  {
     pPlot->setXOption("channel");
     throw;     
  }
  pPlot->xOption(saveXoption);
}

// Additional Declarations