#ifndef DATASET_INCD
#define DATASET_INCD

#include "AAObject.hh"
#include "ana/search/SearchEvent.hh"
#include "astro/Astro.hh"

#include <string>
#include <vector>
#include <algorithm>

using std::string;
using std::vector;

class Model;

struct DataSet : public AAObject
{
  vector<SearchEvent> data;

  // for binned fits: nevents vs logerec vs costhetarec.
  // We actually use the vectors nevents, mean_eproxy, mean_cotheta
  // so we can use the real mean and costheta rather than the bincenter and
  // so that we can skip bins with zero entries fast.

  TH2D* H_nevents; 
  vector<int> nevents; //!
  vector<double> mean_eproxy, mean_costheta; //!


  int id;
  int seed;
  string comment;
  
  DataSet() : H_nevents(nullptr) {}

  string __str__() const;
  
  void prepare_histogram();

  void report_histogram(bool org = false );


  ~DataSet()
  {
    if (  H_nevents  ) delete H_nevents;
    H_nevents = nullptr;
  }

  void clear()
  {
    id = 0;
    seed = 0;
    comment = "";
    data.clear();
    mean_eproxy.clear();
    mean_costheta.clear();
    nevents.clear();
  }

  /*! this is useful in python */
  void copy( DataSet* other ) {
    *this = *other;
  }


  void set_indices() 
  {
    for (unsigned i=0; i < data.size() ; i++ ) data[i].index = i;
  }


  /*! produce a skymap */  
  SkyMap skymap();

  /*! Add this object to a tree */
  void add_to_tree( TTree* T )
  {
    T -> Branch( "data", this, 32000, 4 );
  }

  /*! debug information for model */
  void prnt(Model& M); 


  DataSet select_cone( const EquatorialCoords& center, double conesize ) {

    DataSet r;
    for(auto& e : data ) 
    {
      if (center.distance( e.coordinates ) < conesize ) r.data.push_back( e );
    }

    r.comment = comment + " -> select_cone " + str(conesize ) + " " + center.__str__();
    r.set_indices();
    return r;
  }

  DataSet select_band( double declination, double band_width ) {
    DataSet r;
      for(auto& e : data ) 
        if ( abs( e.coordinates.dec() - declination ) < band_width ) r.data.push_back(e);

    r.comment = comment + " -> select_band" + str(declination) + " " + str( band_width );
    r.set_indices();
    return r;
  }


  vector<DataSet> select_regions( double declination, double conesize ) {

    vector<DataSet> r;

    double d = declination < 0 ? declination - conesize : declination + conesize;


    double l = std::max( 0.0 , 2 * pi * cos ( d ) ) ;

    print ( " select regions ", declination, d, l, conesize );

    int n = l / (2* conesize); // this is how many non-overlappign cookies we can cut at the given decl
    if ( n < 1 ) n = 1;

    for (int i =0; i<n; i++ ) {
      double ra = i * 2 * pi / n;
      r.push_back( select_cone( EquatorialCoords( ra, declination ), conesize ) );
    }

    return r;
  }

  /*! Delete the search-events for which func is true 
      (can be called with python function) */

  void delete_if( bool (*func) ( SearchEvent& evt ) ) {
    data.erase( std::remove_if( data.begin() , data.end(), func ), data.end() );
  }


  ClassDef( DataSet, 1 );
};

namespace cling
{
  inline std::string printValue(const DataSet* t)
  {
    return t->__str__();
  }
}






#endif