#include "DEventProcessor_bcal_mc_timewalk.h"

#include "DANA/DApplication.h"
#include "BCAL/DBCALGeometry.h"
#include "BCAL/DBCALTDCHit.h"
#include "BCAL/DBCALHit.h"
#include "BCAL/DBCALSiPMSpectrum.h"
#include "DHistogram.h"
#include "units.h"

#define DPRINT(x) cout << #x "=" << x << endl

// The executable should define the ROOTfile global variable. It will
// be automatically linked when dlopen is called.
extern TFile *ROOTfile;

// Routine used to create our DEventProcessor
extern "C"{
  void InitPlugin(JApplication *app){
    InitJANAPlugin(app);
    app->AddProcessor(new DEventProcessor_bcal_mc_timewalk());
  }
} // "C"


//------------------
// init
//------------------
jerror_t DEventProcessor_bcal_mc_timewalk::init(void)
{

  for (int i=0; i<4; i++) {
    for (int j=0; j<2; j++) {
      stringstream histo_name;
      histo_name << "signal_" << i << "_" << j;
      signal[i][j] = new TH1D(histo_name.str().c_str(),histo_name.str().c_str(),4000,-100.,300.);
    }
  }

  //set up the tree
  bcal_timewalk_tree = new TTree("bcal_timewalk_tree","bcal_timewalk_tree");

  bcal_timewalk_tree->Branch("E",&m_E,"E/F");
  bcal_timewalk_tree->Branch("first_signal_time",&m_first_signal_time,"first_signal_time/F");
  bcal_timewalk_tree->Branch("last_signal_time",&m_last_signal_time,"last_signal_time/F");
  bcal_timewalk_tree->Branch("t_tdc",&m_t_tdc,"t_tdc/F");
  bcal_timewalk_tree->Branch("t_adc",&m_t_adc,"t_adc/F");
  bcal_timewalk_tree->Branch("module",&m_module,"module/I");
  bcal_timewalk_tree->Branch("layer",&m_layer,"layer/I");
  bcal_timewalk_tree->Branch("sector",&m_sector,"sector/I");
  bcal_timewalk_tree->Branch("end",&m_end,"end/I");

  return NOERROR;
}

//------------------
// brun
//------------------
jerror_t DEventProcessor_bcal_mc_timewalk::brun(JEventLoop *eventLoop, int runnumber)
{
  return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t DEventProcessor_bcal_mc_timewalk::evnt(JEventLoop *loop, int eventnumber)
{
  vector<const DBCALSiPMSpectrum*> spectrums;
  vector<const DBCALTDCHit*> tdc_hits;
  vector<const DBCALHit*> adc_hits;
  loop->Get(spectrums);
  loop->Get(tdc_hits);
  loop->Get(adc_hits);

  map<readout_channel,DHistogram> summed_histo_map;
  map<readout_channel,const DBCALTDCHit*> tdc_hit_map;
  map<readout_channel,const DBCALHit*> adc_hit_map;

  for (int i=0; i < (int)spectrums.size(); i++) {
    int readout_cellId = DBCALGeometry::fADCId(spectrums[i]->module,spectrums[i]->layer,spectrums[i]->sector);
    readout_channel chan(readout_cellId,spectrums[i]->end);

    //Create the histogram if need
    //Since DHistogram default constructor is private, cannot use
    //map::operator[], instead have to use insert() and read()
    if (summed_histo_map.count(chan)==0) {
      summed_histo_map.insert(make_pair(chan,DHistogram(4000,-100.0,300.0)));
    }
    summed_histo_map.find(chan)->second.Add(&(spectrums[i]->spectrum));
  }

  for (int i=0; i<(int)adc_hits.size(); i++) {
    readout_channel chan(adc_hits[i]->cellId, adc_hits[i]->end);

    if (adc_hit_map.count(chan) == 0) {
      adc_hit_map.insert(make_pair(chan,adc_hits[i]));
    } else { //If we have more than one hit in this cell, only take the highest energy one
      if (adc_hits[i]->E > adc_hit_map[chan]->E) {
        adc_hit_map[chan] = adc_hits[i];
      }
    }
  }

  for (int i=0; i<(int)tdc_hits.size(); i++) {
    readout_channel chan(tdc_hits[i]->cellId, tdc_hits[i]->end);

    //if we have more than one hit in the same cell, we will just use whichever
    //one comes first
    if (tdc_hit_map.count(chan) == 0) {
      tdc_hit_map.insert(make_pair(chan,tdc_hits[i]));
    }
  }

  LockState();

  for (map<readout_channel,DHistogram>::const_iterator map_iter = summed_histo_map.begin();
       map_iter != summed_histo_map.end();
       ++map_iter) {
    readout_channel chan = map_iter->first;
    DHistogram summed_spectrum = map_iter->second;

    TH1D *temp_histo = summed_spectrum.MakeTH1D("temp","temp");
    int layer = DBCALGeometry::layer(chan.cellId);
    int end = chan.end == DBCALGeometry::kUpstream ? 0 : 1;
    signal[layer-1][end]->Add(temp_histo);
    delete temp_histo;

    //nothing to do if there's no hit in this cell
    if (adc_hit_map.count(chan)==0) continue;
    const DBCALHit* adc_hit = adc_hit_map[chan];

    const DBCALTDCHit* tdc_hit = NULL;
    if (tdc_hit_map.count(chan)==1) tdc_hit = tdc_hit_map[chan];

    m_E = adc_hit->E;
    m_t_adc = adc_hit->t;
    m_module = adc_hit->module;
    m_layer = adc_hit->layer;
    m_sector = adc_hit->sector;
    m_end = (int)adc_hit->end;

    if(tdc_hit != NULL) {
      m_t_tdc = tdc_hit->t;
    } else {
      m_t_tdc = 0.0;
    }

    m_first_signal_time = summed_spectrum.GetBinCenter(summed_spectrum.FindFirstNonZeroBin());
    m_last_signal_time = summed_spectrum.GetBinCenter(summed_spectrum.FindLastNonZeroBin());

    bcal_timewalk_tree->Fill();
  }

  UnlockState();

  return NOERROR;
}

//------------------
// erun
//------------------
jerror_t DEventProcessor_bcal_mc_timewalk::erun(void)
{
  // Any final calculations on histograms (like dividing them)
  // should be done here. This may get called more than once.
  return NOERROR;
}

//------------------
// fini
//------------------
jerror_t DEventProcessor_bcal_mc_timewalk::fini(void)
{
  return NOERROR;
}