#ifndef Likelihood_C
#define Likelihood_C

#include <TArrow.h>
#include <TPad.h>
#include <TCanvas.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TMath.h>
#include <TFile.h>
#include <TTree.h>
#include <TH1.h>
#include <TObjArray.h>
#include <TVector.h>

//#include <iterator>
#include <vector>
#include <algorithm>
#include <cassert>
#include <iostream>

#include "Geom.C"

namespace Likelihood
{
// double likelihood(const std::vector<int> & clusters, const TH1 * const h)
// { // This function assumes that the clusters are sorted in 
//   // increasing order.  It needs it for efficient iteration.
//   assert(std::is_sorted(clusters.begin(),clusters.end()));
//   int i_c = 0;
//   const int Nbins = h->GetNbinsX();
//   double result = 1; // multiplicative identity
//   for(int i_h = 0; i_h < Nbins; i_h++)
//     { // iterate over all the bins of the histogram
//       // and determine the bin boundaries.
//       const double bincontent = h->GetBinContent(i_h);
//       const double lowx = h->GetBinLowEdge(i_h);
//       const double highx = lowx+h->GetBinWidth(i_h);
//       int inbin = 0;
//       while(lowx < clusters[i_c] and clusters[i_c] <= highx)
//       { // count the multiplicity of the clusters found in this bin
//         inbin++;
//         i_c++;
//       }
//       if(bincontent == 0) continue;

//       // multiply the probability of having that many clusters
//       // in this bin.
//       result *= TMath::Poisson(inbin,bincontent);
//     }
//   return result;
// }

// double LnPoisson(double x, double mu)
// {
//   if(x < 0) return -1*TMath::Infinity();
//   else if (x == 0.0) return -1*mu;
//   else
//   {
//     return x*TMath::Log(mu)-mu-TMath::LnGamma(x+1.0);
//   }
// }

// Less conservative limits:
const std::pair<const double,const double> limits{0.1,0.55};
// More conservative limits:
//const std::pair<const double,const double> limits{0.15,0.53};

double nll2(const std::vector<int> & clusters_v, const TH1* const h,bool quiet = true)
{ // first make an std::map that contains the multiplicities of the
  // clusters in each bin.
  std::map<int,int> clusters;
  double result = 0;
  for(auto c : clusters_v)
  { // std::map::operator[] on an element that doesn't exist
    // automatically creates an element and initializes it, so it
    // starts at 0.
    clusters[c]++;
  }
  // now iterate over the clustertime:multiplicity pairs
  for(const auto &kv : clusters)
  { // kv.first is the cluster time,
    // kv.second is the multiplicity-> ALWAYS GREATER THAN 0 BY CONSTRUCTION.
    const int n_bin = h->FindFixBin(kv.first);
    
    // Old implementation
    //const double likelihood = TMath::Poisson(kv.second,h->GetBinContent(n_bin));
    //if(likelihood == 0) continue;
    //const double loglikelihood = TMath::Log(likelihood);
    //result -= loglikelihood;

    // New implementation
    double mu = h->GetBinContent(n_bin);
    const int x = kv.second;
    assert(x > 0);

    if(mu == 0) //continue; // Poisson(x,0) returns 0, so nll would be -inf.
    { // Don't give up yet, look for adjacent non-zero bins.
      int j = 1;
      int n = 1;
      while(mu == 0)
      {
        mu = (mu*n + h->GetBinContent(n_bin+j) + h->GetBinContent(n_bin-j))/(n+2);
        n += 2;
        j++;
      }
      if(n > h->GetNbinsX()/4)
      {
        std::cout << "Averaged over " << n << " adjacent bins for Likelihood calculation." << std::endl;
      }
    }
    //else std::cout << "No averaging necessary." << std::endl;
    //if(x == 0.0) result += mu; // Avoid the expensive LnGamma call. // USELESS x > 0 always.
    const double quantity = x*TMath::Log(mu)-mu-TMath::LnGamma(x+1.0);
    if(not quiet)
    {
      std::cout << TString::Format("%d,%d,%g,%d,%g",kv.first,n_bin,mu,x,quantity) << std::endl;
    }
    result -= quantity;
    
  }
  return result;
}

double nll(const std::vector<int> & clusters, const TH1 * const h)
{ // This function assumes that the clusters are sorted in 
  // increasing order.  It needs it for efficient iteration.
  assert(std::is_sorted(clusters.begin(),clusters.end()));
  int i_c = 0;
  const int n_clusters = clusters.size();
  const int Nbins = h->GetNbinsX();
  double result = 0; // additive identity
  // Skip clusters arriving before the first bin.
  while(i_c < n_clusters and clusters[i_c] < h->GetBinLowEdge(1)) i_c++;
  for(int i_h = 1; i_h < Nbins+1; i_h++)
  { // iterate over all the bins of the histogram
    // and determine the bin boundaries.
    const double bincontent = h->GetBinContent(i_h);
    const double lowx = h->GetBinLowEdge(i_h);
    const double highx = lowx+h->GetBinWidth(i_h);
    int inbin = 0;
    while(i_c < n_clusters and lowx <= clusters[i_c] and clusters[i_c] < highx)
    { // count the multiplicity of the clusters found in this bin
      inbin++;
      i_c++;
    }
    // We skip empty bins which will return 0 probability
    // for non-zero clusters, and 1 for no clusters. I have to skip after
    // the above loop so that clusters that end up in empty bins are counted
    // properly.
    if(bincontent == 0) continue;

    // add the negative log-likelihood of having that many clusters
    // in this bin.
    const double contribution = TMath::Log(TMath::Poisson(inbin,bincontent));
    result -= contribution;
  }
  return result;
}

double kolmogorov(const std::vector<int> & clusters, const TH1 & h)
{
  if(clusters.size() == 0 or h.GetEntries() == 0) return 0;
  const TAxis * const hxax = h.GetXaxis();
  const double * const hbins = hxax->GetXbins()->GetArray();
  const int n_hbins = hxax->GetNbins();

  TH1D hclu = hbins ? TH1D("hclu","hclu",n_hbins,hbins) : 
                      TH1D("hclu","hclu",n_hbins,hxax->GetXmin(),hxax->GetXmax());
  for(auto & clu : clusters)
  {
    hclu.Fill(clu);
  }
  return h.KolmogorovTest(&hclu,"UO");
}

void test(const std::vector<int> * v = nullptr)
{
  const double bincontents[] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001768972227136034,
  0.0, 0.0007075888908544136, 0.0056607111268353084, 0.006898991685830532, 0.0070758889085441356,
  0.005483813904121705, 0.006014505572262515, 0.005837608349548912, 0.0063683000176897216, 0.006014505572262515,
  0.005130019458694499, 0.007252786131257739, 0.005837608349548912, 0.008491066690252963, 0.006722094463116929,
  0.0070758889085441356, 0.0077834777993985495, 0.0056607111268353084, 0.0063683000176897216, 0.005837608349548912,
  0.005837608349548912, 0.0070758889085441356, 0.0056607111268353084, 0.0070758889085441356, 0.007960375022112153,
  0.004422430567840085, 0.0056607111268353084, 0.0056607111268353084, 0.006898991685830532, 0.007252786131257739,
  0.004599327790553688, 0.005837608349548912, 0.007960375022112153, 0.006722094463116929, 0.005306916681408102,
  0.006014505572262515, 0.005483813904121705, 0.0070758889085441356, 0.005837608349548912, 0.0063683000176897216,
  0.004776225013267291, 0.0070758889085441356, 0.004068636122412878, 0.005130019458694499, 0.006191402794976118,
  0.005483813904121705, 0.006191402794976118, 0.007429683353971342, 0.0063683000176897216, 0.005483813904121705,
  0.005306916681408102, 0.006898991685830532, 0.004776225013267291, 0.005130019458694499, 0.005837608349548912,
  0.006014505572262515, 0.0070758889085441356, 0.0056607111268353084, 0.004599327790553688, 0.0056607111268353084,
  0.005130019458694499, 0.0070758889085441356, 0.004776225013267291, 0.004776225013267291, 0.006545197240403326,
  0.0035379444542720678, 0.0056607111268353084, 0.005306916681408102, 0.006191402794976118, 0.004953122235980895,
  0.004776225013267291, 0.004953122235980895, 0.006014505572262515, 0.006014505572262515, 0.007252786131257739,
  0.004776225013267291, 0.005483813904121705, 0.004599327790553688, 0.007960375022112153, 0.004953122235980895,
  0.004245533345126481, 0.0056607111268353084, 0.007429683353971342, 0.006898991685830532, 0.005130019458694499,
  0.0038917388996992748, 0.004953122235980895, 0.006191402794976118, 0.004776225013267291, 0.004245533345126481,
  0.003714841676985671, 0.004953122235980895, 0.005130019458694499, 0.005483813904121705, 0.006014505572262515,
  0.0063683000176897216, 0.0033610472315584645, 0.004422430567840085, 0.0038917388996992748, 0.005837608349548912,
  0.004776225013267291, 0.007606580576684945, 0.0038917388996992748, 0.004422430567840085, 0.0056607111268353084,
  0.006014505572262515, 0.004953122235980895, 0.006898991685830532, 0.005130019458694499, 0.007252786131257739,
  0.006722094463116929, 0.004599327790553688, 0.004245533345126481, 0.005837608349548912, 0.004953122235980895,
  0.004599327790553688, 0.006191402794976118, 0.004776225013267291, 0.005483813904121705, 0.0056607111268353084,
  0.006722094463116929, 0.0038917388996992748, 0.004068636122412878, 0.004599327790553688, 0.003714841676985671,
  0.0063683000176897216, 0.004245533345126481, 0.006191402794976118, 0.005483813904121705, 0.004776225013267291,
  0.006014505572262515, 0.007252786131257739, 0.0056607111268353084, 0.004953122235980895, 0.004068636122412878,
  0.004599327790553688, 0.005483813904121705, 0.0038917388996992748, 0.004953122235980895, 0.004776225013267291,
  0.004422430567840085, 0.004068636122412878, 0.004599327790553688, 0.004599327790553688, 0.0031841500088448608,
  0.004422430567840085, 0.002653458340704051, 0.0035379444542720678, 0.0035379444542720678, 0.004599327790553688,
  0.0028303555634176542, 0.003714841676985671, 0.002299663895276844, 0.0017689722271360339, 0.0030072527861312575,
  0.0030072527861312575, 0.002653458340704051, 0.002653458340704051, 0.0021227666725632407, 0.0015920750044224304,
  0.002299663895276844, 0.0015920750044224304, 0.0014151777817088271, 0.002653458340704051,
  0.002299663895276844, 0.0021227666725632407, 0.0012382805589952238, 0.0030072527861312575, 0.0014151777817088271,
  0.0017689722271360339, 0.0019458694498496374, 0.0015920750044224304, 0.0012382805589952238, 0.0017689722271360339,
  0.0019458694498496374, 0.002299663895276844, 0.0008844861135680169, 0.002299663895276844, 0.0008844861135680169,
  0.0014151777817088271, 0.0008844861135680169, 0.0021227666725632407, 0.0019458694498496374, 0.0012382805589952238,
  0.0021227666725632407, 0.0014151777817088271, 0.0019458694498496374, 0.0010613833362816203, 0.0019458694498496374,
  0.0019458694498496374, 0.0012382805589952238, 0.0017689722271360339, 0.0015920750044224304, 0.0017689722271360339,
  0.0015920750044224304, 0.0012382805589952238, 0.0017689722271360339, 0.0007075888908544136, 0.0010613833362816203,
  0.0017689722271360339, 0.002653458340704051, 0.0005306916681408102, 0.0014151777817088271, 0.0015920750044224304,
  0.0010613833362816203, 0.0010613833362816203, 0.0015920750044224304, 0.0007075888908544136, 0.0007075888908544136,
  0.0017689722271360339, 0.0014151777817088271, 0.0010613833362816203, 0.0001768972227136034, 0.0007075888908544136,
  0.0010613833362816203, 0.0005306916681408102, 0.0010613833362816203, 0.0008844861135680169, 0.0019458694498496374,
  0.0005306916681408102, 0.0005306916681408102, 0.0012382805589952238, 0.0003537944454272068, 0.0010613833362816203,
  0.0007075888908544136, 0.0005306916681408102, 0.0021227666725632407, 0.0007075888908544136, 0.0005306916681408102,
  0.0003537944454272068, 0.0010613833362816203, 0.0008844861135680169, 0.0008844861135680169, 0.0007075888908544136,
  0.0003537944454272068, 0.0003537944454272068, 0.0007075888908544136, 0.0012382805589952238, 0.0007075888908544136,
  0.0005306916681408102, 0.0010613833362816203, 0.0005306916681408102, 0.0005306916681408102, 0.0008844861135680169,
  0.0003537944454272068, 0.0, 0.0012382805589952238, 0.0001768972227136034, 0.0008844861135680169,
  0.0007075888908544136, 0.0005306916681408102, 0.0005306916681408102, 0.0008844861135680169, 0.0012382805589952238,
  0.0012382805589952238, 0.0003537944454272068, 0.0003537944454272068, 0.0001768972227136034, 0.0012382805589952238,
  0.0007075888908544136, 0.0003537944454272068, 0.0005306916681408102, 0.0005306916681408102, 0.0003537944454272068,
  0.0007075888908544136, 0.0007075888908544136, 0.0010613833362816203, 0.0001768972227136034, 0.0,
  0.0008844861135680169, 0.0001768972227136034, 0.0005306916681408102, 0.0, 0.0001768972227136034,
  0.0001768972227136034, 0.0003537944454272068, 0.0003537944454272068, 0.0005306916681408102, 0.0003537944454272068,
  0.0005306916681408102, 0.0007075888908544136, 0.0001768972227136034, 0.0003537944454272068, 0.0005306916681408102,
  0.0, 0.0008844861135680169, 0.0003537944454272068, 0.0001768972227136034, 0.0008844861135680169,
  0.0005306916681408102, 0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034,
  0.0001768972227136034, 0.0003537944454272068, 0.0003537944454272068, 0.0001768972227136034, 0.0003537944454272068, 0.0005306916681408102,
  0.0, 0.0001768972227136034, 0.0003537944454272068, 0.0008844861135680169, 0.0001768972227136034,
  0.0003537944454272068, 0.0003537944454272068, 0.0003537944454272068, 0.0005306916681408102, 0.0005306916681408102,
  0.0005306916681408102, 0.0001768972227136034, 0.0003537944454272068, 0.0001768972227136034, 0.0001768972227136034,
  0.0001768972227136034, 0.0, 0.0, 0.0, 0.0001768972227136034, 0.0008844861135680169, 0.0001768972227136034,
  0.0001768972227136034, 0.0001768972227136034, 0.0003537944454272068, 0.0001768972227136034, 0.0,
  0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0005306916681408102, 0.0, 0.0005306916681408102,
  0.0003537944454272068, 0.0003537944454272068, 0.0, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0001768972227136034,
  0.0001768972227136034, 0.0001768972227136034, 0.0003537944454272068, 0.0, 0.0001768972227136034, 0.0005306916681408102,
  0.0001768972227136034, 0.0, 0.0, 0.0001768972227136034, 0.0003537944454272068, 0.0001768972227136034,
  0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0001768972227136034,
  0.0, 0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
  0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034, 0.0,
  0.0, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0,
  0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0,
  0.0003537944454272068, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0001768972227136034, 0.0, 0.0, 0.0, 0.0,
  0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034, 0.0001768972227136034, 0.0,
  0.0, 0.0, 0.0, 0.0001768972227136034, 0.0001768972227136034, 0.0, 0.0001768972227136034, 0.0, 0.0,
  0.0012382805589952238, 0.0, 0.0, 0.0, 0.0005306916681408102, 0.0001768972227136034, 0.0, 0.0,
  0.0001768972227136034, 0.0001768972227136034, 0.0, 0.0};
  TH1D * h = new TH1D("h","h",448,128,1024);
  for(int i = 0; i < 448; i++)
  {
    h->SetBinContent(i,bincontents[i]);
  }
  h->Draw();
  const int clusters_arr[] = {50,200,240,300,500,900};
  std::vector<int> clusters;
  for(auto & c : clusters_arr)
  {
    clusters.push_back(c);
  }
  if(v == nullptr) v = &clusters;
  //std::cout << "Likelihood: " << likelihood(*v,h) << std::endl;
  std::cout << "nll2: " << nll2(*v,h) << std::endl;
  std::cout << "nll: " << nll(*v,h) << std::endl;
  std::cout << "kolmogorov: " << kolmogorov(*v,*h) << std::endl;
}

void test(const std::vector<int> & v) {test(&v);}

// int main(int argc, const char * argv[])
// {
//   std::vector<int> v;
//   for(int i = 1; i < argc; i++)
//   {
//     TString tsarg(argv[i]);
//     assert(tsarg.IsDigit() && "Arguments must all be integers, space-separated.");
//     v.push_back(tsarg.Atoi());
//   }
//   test(v);
// }

  class LikelihoodInfo
  {
  public:
    LikelihoodInfo(){};
    LikelihoodInfo(TString fname){ Init(fname); };

    // // I timed, get_nll_slow is about 58% slower than get_nll
    // double get_nll_slow(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    // {
    //   auto nllcompare = [&clusters](const TObject* h1, const TObject* h2)->bool{
    //     return nll(clusters,static_cast<const TH1D*>(h1)) < nll(clusters,static_cast<const TH1D*>(h2)); 
    //   };
    //   auto min_e = std::min_element(std::begin(*(toas.at(side))),
    //                                 std::end(*(toas.at(side))),
    //                                 nllcompare);
    //   return (*ranges.at(side))[toas.at(side)->IndexOf(*min_e)];
    // };

    // double get_nll_rng(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    // {
    //   double min_value = std::numeric_limits<double>::max();
    //   double min_index = -1;
    //   int index = 0;
    //   for(const auto& tobjp : *(toas.at(side)))
    //   {
    //     const double this_value = nll(clusters,static_cast<const TH1D*>(tobjp));
    //     if(this_value < min_value)
    //     {
    //       min_value = this_value;
    //       min_index = index;
    //     }
    //     index++;
    //   }
    //   return (*ranges.at(side))[min_index];
    //};

    double get_nll(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      double min_value = std::numeric_limits<double>::max();
      double min_index = -1;
      int index = 0;
      auto toa = toas.at(side);
      //for(const auto& tobjp : *(toas.at(side)))
      for(const int N = toa->GetEntries(); index < N; index++)
      {
        const double this_value = nll(clusters,static_cast<const TH1D*>(toa->At(index)));
        if(this_value < min_value)
        {
          min_value = this_value;
          min_index = index;
        }
        //index++;
      }
      return (*ranges.at(side))[min_index];
    };

    // double get_nll2_slow(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    // {
    //   auto nll2compare = [&clusters](const TObject* h1, const TObject* h2)->bool{
    //     return nll2(clusters,static_cast<const TH1D*>(h1)) < nll2(clusters,static_cast<const TH1D*>(h2)); 
    //   };
    //   auto min_e = std::min_element(std::begin(*(toas.at(side))),
    //                                 std::end(*(toas.at(side))),
    //                                 nll2compare);
    //   return (*ranges.at(side))[toas.at(side)->IndexOf(*min_e)];
    // };

    double get_nll2(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH,bool quiet=true) const
    {
      double min_value = std::numeric_limits<double>::max();
      double min_x = TMath::QuietNaN();
      //double min_index = -1;
      int index = 0;
      auto toa = toas.at(side);
      auto range = ranges.at(side);
      for(const int N = toa->GetEntries(); index < N; index++)
      {
        // Enforce the upper and lower limits.
        const double curr_x = (*range)[index];
        if(curr_x < limits.first or curr_x > limits.second) continue;

        const double this_value = nll2(clusters,static_cast<const TH1D*>(toa->At(index)),quiet);
        if(this_value < min_value)
        {
          min_value = this_value;
          //min_index = index;
          min_x = curr_x;
        }
      }
      //return (*range)[min_index];
      return min_x;
    };


    // double get_kolmo_slow(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    // {
    //   auto kcompare = [&clusters](const TObject* h1, const TObject* h2)->bool{
    //     return kolmogorov(clusters,static_cast<const TH1D*>(h1)) < kolmogorov(clusters,static_cast<const TH1D*>(h2)); 
    //   };
    //   auto max_e = std::max_element(std::begin(*(toas.at(side))),
    //                                 std::end(*(toas.at(side))),
    //                                 kcompare);
    //   return (*ranges.at(side))[toas.at(side)->IndexOf(*max_e)];
    // };

    double get_kolmo(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      double max_value = 0;
      double max_index = -1;
      int index = 0;
      auto toa = toas.at(side);
      //for(const auto& tobjp : *(toas.at(side)))
      //auto oldErrorIgnoreLevel = gErrorIgnoreLevel;
      //gErrorIgnoreLevel = kError+1;
      for(const int N = toa->GetEntries(); index < N; index++)
      {
        const double this_value = kolmogorov(clusters,*static_cast<const TH1D*>(toa->At(index)));
        if(this_value > max_value)
        {
          max_value = this_value;
          max_index = index;
        }
        //index++;
      }
      //gErrorIngoreLevel = oldErrorIgnoreLevel;
      return (*ranges.at(side))[max_index];
    };


    void Init(TString fname)
    {
      toas.clear();
      ranges.clear();
      TFile f(fname,"READ");
      assert(!f.IsZombie() && "File name to initialize LikelihoodInfo not valid.");
      for(auto & side : {Geom::LEFT,Geom::RIGHT,Geom::BOTH})
      {
        const auto ss = Geom::SideStrings[side];
        const auto toa_name = (TString("toa_")+ss).Data();
        toas[side] = (TObjArray*)f.Get(toa_name);
        const auto range_name = (TString("ranges_")+ss).Data();
        ranges[side] = (TVectorD*)f.Get(range_name);  
      }
    };

    const TObjArray& get_toa(Geom::Side side) const
    {
      return *(toas.at(side));
    };

    const TVectorD& get_range(Geom::Side side) const
    {
      return *(ranges.at(side));
    };

    TGraph * graph_nll(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      TVectorD nlls(toas.at(side)->GetEntries());
      int i = 0;
      auto toa = toas.at(side);
      //for(const auto& op : *(toas.at(side)))
      for(const int N = toa->GetEntries(); i < N; i ++)
      {
        nlls[i] = nll(clusters,static_cast<const TH1D*>(toa->At(i)));
        //i++;
      }
      return new TGraph(*ranges.at(side),nlls);
    };

    TGraph * graph_nll2(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      TVectorD nlls(toas.at(side)->GetEntries());
      int i = 0;
      auto toa = toas.at(side);
      //for(const auto& op : *(toas.at(side)))
      for(const int N = toa->GetEntries(); i < N; i ++)
      {
        nlls[i] = nll2(clusters,static_cast<const TH1D*>(toa->At(i)));
        //i++;
      }
      return new TGraph(*ranges.at(side),nlls);
    };

    TGraph * graph_kolmo(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      TVectorD nlls(toas.at(side)->GetEntries());
      int i = 0;
      auto toa = toas.at(side);
      //for(const auto& op : *(toas.at(side)))
      for(const int N = toa->GetEntries(); i < N; i ++)
      {
        nlls[i] = 1.0-kolmogorov(clusters,*static_cast<const TH1D*>(toa->At(i)));
        //i++;
      }
      return new TGraph(*ranges.at(side),nlls);
    };

    TCanvas * graph_all(const std::vector<int> & clusters, Geom::Side side = Geom::BOTH) const
    {
      TCanvas * c1 = new TCanvas();
      c1->Divide(1,2);
      c1->cd(1);
      TGraph * g_nll = graph_nll(clusters,side);
      TGraph * g_nll2 = graph_nll2(clusters,side);

      g_nll2->SetLineColor(kRed);
      
      TMultiGraph * mg = new TMultiGraph("mg","Comparing Likelihoods");
      mg->Add(g_nll);
      mg->Add(g_nll2);
      
      mg->Draw("AL");
      mg->GetYaxis()->SetTitle("Negative Log Likelihood");
      mg->GetXaxis()->SetTitle("Track distance (cm)");
      

      TArrow * nllline = new TArrow();
      auto x0_nll = get_nll(clusters,side);
      //auto mg_x = g_nll->Eval(best_nll);
      auto g_min = g_nll->GetYaxis()->GetXmin();
      auto g_max = g_nll->GetYaxis()->GetXmax();
      //std::cout << x0_nll << ", " << g_min << ", " << g_max << std::endl;
      nllline->DrawArrow(x0_nll,g_min,x0_nll,g_min + (g_max-g_min)*0.1,0.01);

      TArrow * nll2line = new TArrow();
      nll2line->SetLineColor(kRed);
      nll2line->SetFillColor(kRed);
      auto x0_nll2 = get_nll2(clusters,side);
      auto g2_min = g_nll2->GetYaxis()->GetXmin();
      auto g2_max = g_nll2->GetYaxis()->GetXmax();
      nll2line->DrawArrow(x0_nll2,g2_min,x0_nll2,g2_min + (g2_max-g2_min)*0.1,0.01);

      TLine * limitline = new TLine();
      limitline->DrawLine(limits.first,g2_min,limits.first,g2_max);
      limitline->DrawLine(limits.second,g2_min,limits.second,g2_max);

      gPad->Modified();
      gPad->Update();

      c1->cd(2);
      TGraph * g_kolmo = graph_kolmo(clusters,side);
      g_kolmo->Draw("AL");
      g_kolmo->GetYaxis()->SetTitle("KS Test Statistic");
      g_kolmo->GetXaxis()->SetTitle("Track distance (cm)");
      //gPad->SetLogy();
      TArrow * kolmoline = new TArrow();
      auto x0_kolmo = get_kolmo(clusters,side);
      auto gk_min = g_kolmo->GetYaxis()->GetXmin();
      auto gk_max = g_kolmo->GetYaxis()->GetXmax();
      kolmoline->DrawArrow(x0_kolmo,gk_min,x0_kolmo,gk_min + (gk_max-gk_min)*0.1,0.01);

      gPad->Modified();
      gPad->Update();

      // c1->cd(3);
      // auto b
      // TH1D * hclu = new TH1D("hclu","hclu",)
      // toas.at(side)->At()
      g_nll->SetBit(kCanDelete);
      g_nll2->SetBit(kCanDelete);
      g_kolmo->SetBit(kCanDelete);
      mg->SetBit(kCanDelete);
      nllline->SetBit(kCanDelete);
      nll2line->SetBit(kCanDelete);
      kolmoline->SetBit(kCanDelete);
      limitline->SetBit(kCanDelete);
      return c1;
    }    

  //private:
    std::map<Geom::Side,TObjArray*> toas; // TObjArrays of TH1D* with cluster distros.
    std::map<Geom::Side,TVectorD *> ranges;

  };

  void FourWayGraph(const LikelihoodInfo & lhi, 
                    const std::vector<int> & clusters,
                    double trad_dist,
                    Geom::Side theside = Geom::BOTH,
                    TString savename = "outfiles/plots/fourway.pdf")
  {

    // Graph of nll values for the clusters.
    std::unique_ptr<TGraph> nll2_graph(lhi.graph_nll2(clusters,theside));
    nll2_graph->GetXaxis()->SetTitle("Distance from wire (cm)");
    nll2_graph->GetYaxis()->SetTitle("Negative Log Likelihood");
    nll2_graph->SetTitle("Red arrow = traditional tracking, black arrow = minimum NLL");

    // Now get he best nll distance, and extract the TH1D* for that distance.
    const double nll_dist = lhi.get_nll2(clusters,theside);
    const TVectorD * const ranges_tv = lhi.ranges.at(theside);
    const double * const ranges_arr = ranges_tv->GetMatrixArray();
    const int ranges_N = ranges_tv->GetNoElements();
    const int best_idx = std::upper_bound(ranges_arr,ranges_arr+ranges_N,nll_dist) - 1 - ranges_arr;
    const std::unique_ptr<TH1D> best_hist(static_cast<TH1D*>(lhi.toas.at(theside)->At(best_idx)->Clone("best_hist")));
    const double check1 = nll2(clusters,best_hist.get()); 
    //const double check2 = *std::min_element(nll2_graph->GetY(),nll2_graph->GetY()+nll2_graph->GetN());
    //assert(check1 == check2);
    //best_hist->SetTitle("Maximum Likelihood Cluster Distribution");

    // Now do the same for the TH1D* for the distance found by traditional tracking.
    const int trad_idx = std::upper_bound(ranges_arr,ranges_arr+ranges_N,trad_dist) - 1 - ranges_arr;
    //TH1D * trad_hist = nullptr;
    std::unique_ptr<TH1D> trad_hist(nullptr);
    if(trad_idx < ranges_N-1)
    {
      trad_hist.reset(static_cast<TH1D*>(lhi.toas.at(theside)
                      ->At(trad_idx)->Clone("trad_hist")));    
    }
    //if(trad_hist) trad_hist->SetTitle("Traditional Tracking Cluster Distribution");

    // Now construct a "histogram" of the cluster times, just like we used to do
    // for the Kolmogorov Test.
    const TAxis * const hxax = best_hist->GetXaxis();
    const double * const hbins = hxax->GetXbins()->GetArray();
    const int n_hbins = hxax->GetNbins();
    TString htitle("Actual Clusters Found");
    TH1D hclu = hbins ? TH1D("hclu",htitle,n_hbins,hbins) : 
                        TH1D("hclu",htitle,n_hbins,hxax->GetXmin(),hxax->GetXmax());
    for(auto & clu : clusters)
    {
      hclu.Fill(clu);
    }

    TLine line;
    TArrow arrow;
    auto g_min = nll2_graph->GetYaxis()->GetXmin();
    auto g_max = nll2_graph->GetYaxis()->GetXmax();

    // Draw stuff.
    TCanvas c1;
    c1.Divide(2,2);
    c1.cd(1);
    hclu.Draw();

    c1.cd(2);
    nll2_graph->Draw("AL");
    arrow.DrawArrow(nll_dist,g_min,nll_dist,g_min + (g_max-g_min)*0.1,0.01);
    arrow.SetLineColor(kRed);
    arrow.SetFillColor(kRed);
    arrow.DrawArrow(trad_dist,g_min,trad_dist,g_min + (g_max-g_min)*0.1,0.01);

    line.DrawLine(limits.first,g_min,limits.first,g_max);
    line.DrawLine(limits.second,g_min,limits.second,g_max);

    c1.cd(3);
    best_hist->Draw();

    c1.cd(4);
    if(trad_hist) trad_hist->Draw();


    c1.SaveAs(savename);
  }

}

#endif