#include "Det.hh"
#include "TH1D.h"

// This example contains c++ code to fill histograms with a delta-distribution
// from timeslice data, for all pars of pmts on all doms.
// The code uses pure c-style arrays for speed.


inline void init_histograms( vector<TH1D*>& HH, int domid  )
{
  HH.resize(31*31);
  
  for (int c1=0; c1<31; ++c1)
      for (int c2=c1+1; c2<31;c2++)
	{
	  const int id = c1*31+c2;
	  HH[id] = new TH1D(Form("delta_t_%d_%d+%d",domid, c1, c2 ),"h",51,-25,25);
	}
}

// Process hits, increasing hit_index untill dom_id changes.
//
// Because of how they are read from the timeslice, the hits for one dom
// are guarranteed to be continuous. Furthermore, for each given PMT, the
// hits are time-sorted. Both properties are exploited.
//
// The detector objects should hold the rates for each Pmt.

inline void proc_one_dom( vector<Hit>& hits, unsigned& hit_index, TH1D** H , Det& det )
{
  static double M[31][1000];
  int N[31] = {0}; 

  const int domid = hits[hit_index].dom_id;

  // We keep only pmts that have between 10 and 10000 (L0) hits according
  // to the summary data. Also store these rates in an array to save map-lookups.
  static int rates[31];
  static bool channel_ok[31];
  
  for (int i=0; i<31; ++i)
    {
      rates[i]      = det.doms[domid].pmts[i].rate;
      channel_ok[i] = !( rates[i]<10 or rates[i]>10000 );
    }

  static const double f  = 0.10 * H[1]->GetBinWidth(0)*1e-9;
  static const int overflow = H[1]->GetNbinsX()+1;

  // Hit loop, fill array of hit.t for each PMT.
  for( ;; ++hit_index )
    {
      Hit& h = hits[hit_index];
      if ( h.dom_id != domid ) break; 
      if ( h.t == 0          ) continue ; // what is causing hits with t==0 ??
      
      const unsigned c = h.channel_id; 
      if ( !channel_ok[c] ) continue;

      M[c][ N[c] ] = h.t;

      if ( ++N[c] > 1000 )
	{
	  cout << "too many hits on dom" << domid << " pmt " << c << " dropping remaining hits " << endl;;
	  break;
	}
    }

  // Loop over all combinations of PMTs
  for (int c1=0; c1<31; ++c1)
    {
      if (!channel_ok[c1]) continue;
      
      const int N1 = N[c1]+1;
      const int Y  = c1*31;
      
      for (int c2=c1+1; c2<31;++c2)
	{
	  if (!channel_ok[c2]) continue;
	  
	  const int N2 = N[c2]+1;
	  const int id = Y + c2;
	  
	  // We (ab)use the underflow bin to store the total expected bg
	  // coincidence rate per bin. The overflow bin holds the number of
	  // timeslices where both PMTs were active
	  const double rate = f*rates[c1]*rates[c2];
	  H[id]->AddBinContent(0,        rate   );
	  H[id]->AddBinContent(overflow, 1 );

	  // Loop over all pairs on these two PMTs, using the fact that
	  // the hits are sorted on each PMT.
	  
	  for(int i1(0), i2(0);;)
	    {
	      const double dt =  M[c2][i2] - M[c1][i1];
	      if      ( dt > 20 ) {if (++i1==N1) break;}
	      else if ( dt <-20 ) {if (++i2==N2) break;}
	      else
		{
		  H[id]->Fill( dt );
		  if (++i1==N1 or ++i2==N2) break; // assumes every hits contributes only to one L1/pair
		}
	    }	   		  
	}
    } 
}


inline map<int, vector< TH1D* > >&  process_coincidences( vector<Hit>& hits, Det& det )
{ 
  static map<int, vector< TH1D* > > histograms;
  
  hits.push_back( Hit() ); // sentinel with dom_id =0;
 
  unsigned idx=0;
  while (idx < hits.size()-1) 
    {
      int domid = hits[idx].dom_id;
      vector<TH1D*>& HH = histograms[ domid ];
      if ( HH.size() ==0 ) init_histograms( HH, domid );
      proc_one_dom( hits, idx , &(HH[0]), det );
    }
  hits.pop_back(); // remove the sentinel hit
  return histograms;
}


// The pragma's below tell root to generate dictionary information for these types,
// wich is needed to analyse the results of the process_coincidences function

#pragma link C++ class map<int, vector< TH1D* > >;
#pragma link C++ class pair<int, vector< TH1D* > >;
#pragma link C++ class map<int, vector< TH1D* > >::iterator;
#pragma link C++ class vector< TH1D* >;
#pragma link C++ class vector< TH1D* >::iterator;