// $Id$
//
// File: JEventProcessor_timing.cc
// Created: Fri Nov  9 11:58:09 EST 2012
// Creator: wolin (on Linux stan.jlab.org 2.6.32-279.11.1.el6.x86_64 x86_64)
// Edited v1: Eric Pooser 10/09/2014, epoos001@fiu.edu
// Edited v2: David Lawrence 10/11/2014, davidl@jlab.org
// Edited v3: Eric Pooser 10/13/2014, epoos001@fiu.edu
// Edited v4: Eric Pooser, Werner Boeglin 10/17/14, epoos001@fiu.edu, boeglinw@fiu.edu
// Edited v5: Eric Pooser, Werner Boeglin 10/17/14, epoos001@fiu.edu, boeglinw@fiu.edu

#include <stdint.h>
#include <vector>
#include <TMath.h>
#include "JEventProcessor_timing.h"
#include <JANA/JApplication.h>
using namespace std;
using namespace jana;

#include <TDirectory.h>
#include <TH1.h>
#include <TH2.h>
#include <TGraphErrors.h>

#include "START_COUNTER/DSCHit.h"
#include "START_COUNTER/DSCDigiHit.h"
#include "START_COUNTER/DSCTDCDigiHit.h"

#include "BCAL/DBCALDigiHit.h"
#include "BCAL/DBCALTDCDigiHit.h"
#include "BCAL/DBCALHit.h"
#include "BCAL/DBCALTDCHit.h"
#include "BCAL/DBCALPoint.h"
#include "BCAL/DBCALCluster.h"
#include "BCAL/DBCALShower.h"
#include "BCAL/DBCALGeometry.h"

#include "FCAL/DFCALHit.h"
#include "FCAL/DFCALDigiHit.h"
#include "FCAL/DFCALGeometry.h"
#include "FCAL/DFCALCluster.h"
#include "FCAL/DFCALShower.h"

#include "CDC/DCDCHit.h"
#include "CDC/DCDCDigiHit.h"

#include "FDC/DFDCHit.h"
#include "FDC/DFDCWireDigiHit.h"
#include "FDC/DFDCCathodeDigiHit.h"

#include "TOF/DTOFHit.h"
#include "TOF/DTOFDigiHit.h"
#include "TOF/DTOFTDCDigiHit.h"
#include "TTAB/DTranslationTable.h"

#include "PAIR_SPECTROMETER/DPSCDigiHit.h"
#include "PAIR_SPECTROMETER/DPSCTDCDigiHit.h"
#include "PAIR_SPECTROMETER/DPSDigiHit.h"

#include "TAGGER/DTAGHDigiHit.h"
#include "TAGGER/DTAGHTDCDigiHit.h"

#include "TAGGER/DTAGMDigiHit.h"
#include "TAGGER/DTAGMTDCDigiHit.h"

#include "DAQ/DCODAEventInfo.h"
#include "DAQ/DCODAROCInfo.h"
#include "DAQ/DF1TDCHit.h"

// Define some constants
const uint32_t  NCHANNELS     = 30;              // number of scintillator paddles
const float_t   TDC_RES       = 0.057;           // f1TDC resolution (ns)
const float_t   ADC_PT_RES    = 0.0625;          // fADC250 pulse time resolution (ns)
// Define some histogram limits
const uint32_t  PI_MIN         = 0;              // Lower limit of pulse integral histogram
const uint32_t  PI_MAX         = 70000;          // Upper limit of pulse integral histogram  
const uint32_t  PI_BINS        = 700;            // Number of bins in pulse integral histogram
const uint32_t  PED_MIN        = 0;              // Lower limit of pedestal histogram
const uint32_t  PED_MAX        = 400;            // Upper limit of pedestal histogram
const uint32_t  PED_BINS       = 400;            // Number of bins in pedestal histogram
const uint32_t  PT_MIN         = 0;              // Lower limit of pulse time histogram (ns)
const uint32_t  PT_MAX         = 400;            // Upper limit of pulse time histogram (ns)
const uint32_t  PT_BINS        = 400;            // Number of bins in pulse time histogram
const uint32_t  TDC_DHIT_MIN   = 0;              // Lower limit of tdc digihit time histogram (ns)
const uint32_t  TDC_DHIT_MAX   = 4000;           // Upper limit of tdc digihit time histogram (ns)
const uint32_t  TDC_DHIT_BINS  = 400;            // Number of bins in tdc digihit time histogram 
const uint32_t  BCAL_TDC_DHIT_MIN   = 0;         // Lower limit of tdc digihit time histogram (ns)
const uint32_t  BCAL_TDC_DHIT_MAX   = 50;        // Upper limit of tdc digihit time histogram (ns)
const uint32_t  BCAL_TDC_DHIT_BINS  = 100;       // Number of bins in tdc digihit time histogram 
const uint32_t  ADC_MULTI_MIN  = 0;              // Lower limit of adc multiplicity histogram
const uint32_t  ADC_MULTI_MAX  = 10;             // Upper limit of adc multiplicity histogram
const uint32_t  ADC_MULTI_BINS = 10;             // Number of bins in adc multiplicity histogram
const uint32_t  TDC_MULTI_MIN  = 0;              // Lower limit of tdc multiplicity histogram
const uint32_t  TDC_MULTI_MAX  = 10;             // Upper limit of tdc multiplicity histogram
const uint32_t  TDC_MULTI_BINS = 10;             // Number of bins in tdc multiplicity histogram
// Define some variables
int32_t tdiff_int;                               // integer value of the time difference between two hits
float_t tdiff;                                   // float value of the time difference between two hits
int32_t adc_multiplicity_counter[NCHANNELS];     //  array for storing the fADC250 multiplicity of the scintillators
int32_t tdc_multiplicity_counter[NCHANNELS];     //  array for storing the f1TDC multiplicity of the scintillators

// DSCDigiHit 1D histosgram pointers
static TH1I *st_sec_adc_dhit;
static TH1I *st_pi_dhit;
static TH1I *st_pt_dhit;
static TH1I *st_ped_dhit;
static TH1I *st_qf_dhit;
static TH1I *st_nsi_dhit;
static TH1I *st_nsp_dhit;
// DSCHit 1D histosgram pointers
static TH1I *st_sec_hit;
//static TH1I *st_de_hit;
//static TH1I *st_t_hit;
//static TH1I *st_sig_t_hit;
// DSCTDCDigiHit 1D histogram pointers
static TH1I *st_sec_tdc_dhit;
static TH1I *st_tdc_dhit;
// DSCDigiHit 2D histosgram pointers
static TH2I *st_pi_dhit_2d;
static TH2I *st_ped_dhit_2d;
static TH2I *st_pt_dhit_2d;
static TH2I *st_pi_pt_dhit_2d;
static TH2I *st_pi_nsi_dhit_2d;
// DSCHit 2D histosgram pointers
//static TH2I *st_t_hit_2d;
//static TH2I *st_sig_t_hit_2d;
//static TH2I *st_de_hit_2d;
//static TH2I *st_de_t_hit_2d;
// DSCTDCDigiHit 2D histosgram pointers
static TH2I *st_tdc_dhit_2d;
static TH2I *st_pi_tdc_dhit_2d;
static TH2I *st_pt_tdc_dhit_2d;
// Dynamic Array of 1D histogram pointers
TH1I** st_pi_dhit_chan   = new TH1I*[NCHANNELS];
TH1I** st_ped_dhit_chan  = new TH1I*[NCHANNELS];
TH1I** st_tdc_dhit_chan  = new TH1I*[NCHANNELS];
//TH1I** st_t_hit_chan     = new TH1I*[NCHANNELS];
//TH1I** st_sig_t_hit_chan = new TH1I*[NCHANNELS];
TH1I** st_adc_multi_chan = new TH1I*[NCHANNELS];
TH1I** st_tdc_multi_chan = new TH1I*[NCHANNELS];
// Number of events histogram pointer
static TH1I *st_num_events;
// Time difference histograms
static TH1I *st_tdc_tdiff;
static TH1I *st_tdc_tdiff_all;
// Multiplicity histograms
static TH1I *st_adc_multi;
static TH1I *st_tdc_multi;
static TH2I *st_adc_multi_2d;
static TH2I *st_tdc_multi_2d;
static TH2I *st_adc_tdc_multi_2d;

// ST & BCAL histograms
static TH1I *st_bcal_tdc;
static TH1I *st_bcal_tdc_us;
static TH1I *st_bcal_tdc_ds;
static TH1I *st_bcal_tdc_avg;
static TH2I *st_bcal_tdc_2d;
static TH2I *st_bcal_tdc_us_2d;
static TH2I *st_bcal_tdc_ds_2d;
static TH2I *st_bcal_tdc_avg_2d;

// DCODAEventInfo
static TH1I *hnEventInfo;

// DCODAROCInfo
static TH1I *hnROCInfo;
static TH1I *hnROCId;

//////////////////////// KM hists /////////////////////////////////////

/// number of digihits for each detector ///
static TH1I *hnST;
static TH1I *hnBCAL;
static TH1I *hnFCAL;
static TH1I *hnCDC;
static TH1I *hnFDCcathode;
static TH1I *hnFDCanode;
static TH1I *hnTOF;

static TH1I *hnST_TDC;
static TH1I *hnBCAL_TDC;

static TH1I *hnTOF_TDC;

/// correlation of digihits ///

// ST vs BCAL
static TH2I *hADC_ST_BCAL;
static TH2I *hADC_ST_BCAL_up;
static TH2I *hADC_ST_BCAL_down;

static TH2I *hTDC_ST_BCAL;
static TH2I *hTDC_ST_BCAL_up;
static TH2I *hTDC_ST_BCAL_down;

// ST vs FCAL
static TH2I *hADC_ST_FCAL;

// ST vs TOF
static TH2I *hADC_ST_TOF;
static TH2I *hTDC_ST_TOF;

// ST ADC vs ST TDC
static TH2I *hST_ADC_TDC;

// DCodaEventInfo avg time
static TGraph *gEventInfoAvgTime;
static TGraph *gtimestamp[96];
//----------------------------------------------------------------------------------
// Routine used to create our JEventProcessor
extern "C"{
  void InitPlugin(JApplication *app){
    InitJANAPlugin(app);
    app->AddProcessor(new JEventProcessor_timing());
  }
}
//----------------------------------------------------------------------------------
JEventProcessor_timing::JEventProcessor_timing() {
}
//----------------------------------------------------------------------------------
JEventProcessor_timing::~JEventProcessor_timing() {
}
//----------------------------------------------------------------------------------

jerror_t JEventProcessor_timing::init(void) {
  
  japp->RootWriteLock(); //ACQUIRE ROOT LOCK!!

  // Set flag of whether to get raw pulse data
  getRawData = false;
  cout << "Flag for getting/saving raw pulse data is " << (getRawData==true ? "true" : "false") << endl;

  // Create root folder for ST and cd to it, store main dir
  TDirectory *main = gDirectory;
  gDirectory->mkdir("timing")->cd();
  // Book DSCDigiHit 1D histosgrams 
  st_sec_adc_dhit = new TH1I("st_sec_adc_dhit", "ST fADC250 DigiHit Occupancy; Channel Number; fADC250 DigiHits", NCHANNELS, 0.5, NCHANNELS + 0.5);
  st_pi_dhit      = new TH1I("st_pi_dhit", "ST fADC250 Pulse Integral; fADC250 Pulse Integral; Counts ",  PI_BINS, PI_MIN, PI_MAX);
  st_pt_dhit      = new TH1I("st_pt_dhit", "ST fADC250 Pulse Time; fADC250 Pulse Time (ns); Counts", PT_BINS, PT_MIN, PT_MAX);
  st_ped_dhit     = new TH1I("st_ped_dhit", "ST fADC250 Pulse Pedestal; fADC250 Pulse Integral; Counts", PED_BINS, PED_MIN, PED_MAX);
  st_qf_dhit      = new TH1I("st_qf_dhit", "ST Quality Factor; fADC250 Quality Factor; Counts", 4, 0, 4);
  st_nsi_dhit     = new TH1I("st_nsi_dhit", "ST Number of Samples Used in Integral; Number of Samples used in Pulse Integral; Counts", 150, 0, 150);
  st_nsp_dhit     = new TH1I("st_nsp_dhit", "ST Number of Samples Used in Pedestal; Number of Samples used in Pulse Pedestal; Counts", 800,0, 800);
  // Book DSCHit 1D histosgrams
  st_sec_hit   = new TH1I("st_sec_hit", "ST Hit Occupancy; Channel Number; Hits", NCHANNELS, 0.5, NCHANNELS + 0.5);
  //st_de_hit    = new TH1I("st_de_hit", "ST Energy Loss; dE/dx (MeV); Counts", 3000, -10.0, 20.0);
  //st_t_hit     = new TH1I("st_t_hit", "ST f1TDC Hit Time; Hit Time (ns); Counts", 6000, -20.0, 40.0);
  //st_sig_t_hit = new TH1I("st_sig_t_hit", "ST #sigma_{t} Hit Time; Hit Time Uncertainty (ns); Counts", 1000, 0.0, 10.0);
  // Book DSCTDCDigiHit 1D histogram
  st_sec_tdc_dhit = new TH1I("st_sec_tdc_dhit", "ST f1TDC DigiHit Occupancy; Channel Number; f1TDC DigiHits", NCHANNELS, 0.5, NCHANNELS + 0.5);
  st_tdc_dhit     = new TH1I("st_tdc_dhit", "ST f1TDC DigiHit Time; f1TDC DigiHit Time (ns); Counts", TDC_DHIT_BINS, TDC_DHIT_MIN, TDC_DHIT_MAX);
  // Book 2D histograms
  st_pi_dhit_2d     = new TH2I("st_pi_dhit_2d", "ST fADC250 Pulse Integral; Channel Number; fADC250 Pulse Integral", NCHANNELS, 0.5, NCHANNELS + 0.5, PI_BINS, PI_MIN, PI_MAX);
  st_ped_dhit_2d    = new TH2I("st_ped_dhit_2d", "ST fADC250 Pulse Pedestal; Channel Number; fADC250 Pulse Pedestal", NCHANNELS, 0.5, NCHANNELS + 0.5, PED_BINS, PED_MIN, PED_MAX);
  st_pt_dhit_2d     = new TH2I("st_pt_dhit_2d", "ST fADC250 Pulse Time; Channel Number; fADC250 Pulse Time (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, PT_BINS, PT_MIN, PT_MAX);
  //st_t_hit_2d       = new TH2I("st_t_hit_2d", "ST TDC Hit Time; Channel Number; Hit Time (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, 6000, -20.0, 40.0);
  //st_sig_t_hit_2d   = new TH2I("st_sig_t_hit_2d", "ST #sigma_{t} Hit Time; Channel Number; Hit Time Uncertainty (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, 1000, 0.0, 10.0);
  //st_de_hit_2d      = new TH2I("st_de_hit_2d", "ST dE/dx; Channel Number; dE/dx (MeV)", NCHANNELS, 0.5, NCHANNELS + 0.5, 3000, -10.0, 20.0);
  st_tdc_dhit_2d    = new TH2I("st_tdc_dhit_2d", "ST f1TDC DigiHit Time; Channel Number; f1TDC DigiHit Time (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, TDC_DHIT_BINS, TDC_DHIT_MIN, TDC_DHIT_MAX);
  st_pi_pt_dhit_2d  = new TH2I("st_pi_pt_dhit_2d", "ST Pulse Integral & Pulse Time Correlation; fADC250 Pulse Integral; fADC250 Pulse Time (ns);", PI_BINS, PI_MIN, PI_MAX, PT_BINS, PT_MIN, PT_MAX);
  st_pi_nsi_dhit_2d = new TH2I("st_pi_nsi_dhit_2d", "ST Pulse Integral & Number of Samples Correlation; fADC250 Number of Samples Used in Pulse Integral; fADC250 Pulse Integral", 150, 0, 150, PI_BINS, PI_MIN, PI_MAX);
  st_pi_tdc_dhit_2d = new TH2I("st_pi_tdc_dhit_2d", "ST fADC250 Pulse Integral & f1TDC Time DigiHit Correlation; f1TDC DigiHit Time (ns); fADC250 Pulse Integral", TDC_DHIT_BINS, TDC_DHIT_MIN, TDC_DHIT_MAX, PI_BINS, PI_MIN, PI_MAX);
  //st_de_t_hit_2d    = new TH2I("st_de_t_hit_2d", "ST dE/dx & TDC Time Hit Correlation; Hit Time (ns); dE/dx (MeV)", 6000, -20., 40., 3000, -10., 20.);
  st_pt_tdc_dhit_2d = new TH2I("st_pt_tdc_dhit_2d", "ST Pulse Time & f1TDC Time DigiHit Correlation; f1TDC DigiHit Time (ns); fADC250 Pulse Time (ns)", TDC_DHIT_BINS, TDC_DHIT_MIN, TDC_DHIT_MAX, PT_BINS, PT_MIN, PT_MAX);
  // Create directoy for the individual paddle histograms
  gDirectory->mkdir("st_channel_histos")->cd();
  // Book dynamic array of 1D histograms
  for(unsigned int i = 0; i < NCHANNELS; i++)
    {
      st_pi_dhit_chan[i]   = new TH1I(Form("st_pi_dhit_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; fADC250 Pulse Integral; Counts", i+1, 0+12*i, 12+12*i), PI_BINS, PI_MIN, PI_MAX);
      st_ped_dhit_chan[i]  = new TH1I(Form("st_ped_dhit_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; fADC250 Pulse Pedestal; Counts", i+1, 0+12*i, 12+12*i), PED_BINS, PED_MIN, PED_MAX);
      st_tdc_dhit_chan[i]  = new TH1I(Form("st_tdc_dhit_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; f1TDC DigiHit Time (ns); Counts", i+1, 0+12*i, 12+12*i), TDC_DHIT_BINS, TDC_DHIT_MIN, TDC_DHIT_MAX );
      //st_t_hit_chan[i]     = new TH1I(Form("st_t_hit_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; Hit Time (ns); Counts", i+1, 0+12*i, 12+12*i), 12000, -20, 100);
      //st_sig_t_hit_chan[i] = new TH1I(Form("st_sig_t_hit_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; Hit Time Uncertainty (ns); Counts", i+1, 0+12*i, 12+12*i), 1000, 0.0, 10.0);
      st_adc_multi_chan[i] = new TH1I(Form("st_adc_multi_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; fADC250 Multiplicty; Counts", i+1, 0+12*i, 12+12*i), ADC_MULTI_BINS, ADC_MULTI_MIN + 0.5, ADC_MULTI_MAX + 0.5);
      st_tdc_multi_chan[i] = new TH1I(Form("st_tdc_multi_chan_%i", i+1), Form("Channel %i, #phi #in [%i^{#circ}, %i^{#circ}]; f1TDC Multiplicty; Counts", i+1, 0+12*i, 12+12*i), TDC_MULTI_BINS, TDC_MULTI_MIN + 0.5, TDC_MULTI_MAX + 0.5);
    }
  gDirectory->cd("../");

  // Book time difference histograms
  st_tdc_tdiff     = new TH1I("st_tdc_tdiff", "ST f1TDC Time Difference; Time Difference (ns); Counts", 400, -10.0, 10.0);
  st_tdc_tdiff_all = new TH1I("st_tdc_tdiff_all", "ST f1TDC Time Difference; Time Difference (ns); Counts", 1800, -300.0, 300.0);
  // Book multiplicity histograms
  st_tdc_multi        = new TH1I("st_tdc_multi", "ST f1TDC Multiplicity; Number of ST f1TDC Hits; Counts", TDC_MULTI_BINS, TDC_MULTI_MIN + 0.5, TDC_MULTI_MAX + 0.5);
  st_adc_multi        = new TH1I("st_adc_multi", "ST fADC250 Multiplicity; Number of ST fADC250 Hits; Counts", ADC_MULTI_BINS, ADC_MULTI_MIN + 0.5, ADC_MULTI_MAX + 0.5);
  st_adc_multi_2d     = new TH2I("st_adc_multi_2d", "ST fADC250 Multiplicity; Channel Number; fADC250 Multiplicty", NCHANNELS, 0.5, NCHANNELS + 0.5, ADC_MULTI_BINS, ADC_MULTI_MIN + 0.5, ADC_MULTI_MAX + 0.5);
  st_tdc_multi_2d     = new TH2I("st_tdc_multi_2d", "ST f1TDC Multiplicity; Channel Number; f1TDC Multiplicity", NCHANNELS, 0.5, NCHANNELS + 0.5, TDC_MULTI_BINS, TDC_MULTI_MIN + 0.5, TDC_MULTI_MAX + 0.5);
  st_adc_tdc_multi_2d = new TH2I("st_adc_tdc_multi_2d", "ST fADC250 vs. f1TDC Multiplicity; f1TDC Multiplicity; fADC250 Multiplicity", TDC_MULTI_BINS, TDC_MULTI_MIN + 0.5, TDC_MULTI_MAX + 0.5, ADC_MULTI_BINS, ADC_MULTI_MIN + 0.5, ADC_MULTI_MAX + 0.5);

  // Book BCAL and ST histograms
  st_bcal_tdc        = new TH1I("st_bcal_tdc", "Time Difference of ST & BCAL TDC DigiHit Times; BCAL_{TDC} - ST_{TDC} (ns); Counts", BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_us     = new TH1I("st_bcal_tdc_us", "Time Difference of ST & Upstream BCAL TDC DigiHit Times; BCAL_{TDC} - ST_{TDC} (ns); Counts", BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_ds     = new TH1I("st_bcal_tdc_ds", "Time Difference of ST & Downstream BCAL TDC DigiHit Times; BCAL_{TDC} - ST_{TDC} (ns); Counts", BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_avg    = new TH1I("st_bcal_tdc_avg", "Time Difference of ST & Average BCAL TDC DigiHit Times; BCAL_{TDC} - ST_{TDC} (ns); Counts", BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_2d     = new TH2I("st_bcal_tdc_2d", "Time Difference of ST & BCAL TDC DigiHit Times; Channel Number; BCAL_{TDC} - ST_{TDC} (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_us_2d  = new TH2I("st_bcal_tdc_us_2d", "Time Difference of ST & BCAL Upstream TDC DigiHit Times; Channel Number; BCAL_{TDC} - ST_{TDC} (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_ds_2d  = new TH2I("st_bcal_tdc_ds_2d", "Time Difference of ST & BCAL Downstream TDC DigiHit Times; Channel Number; BCAL_{TDC} - ST_{TDC} (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);
  st_bcal_tdc_avg_2d = new TH2I("st_bcal_tdc_avg_2d", "Time Difference of ST & BCAL Average TDC DigiHit Times; Channel Number; BCAL_{TDC} - ST_{TDC} (ns)", NCHANNELS, 0.5, NCHANNELS + 0.5, BCAL_TDC_DHIT_BINS, BCAL_TDC_DHIT_MIN, BCAL_TDC_DHIT_MAX);

  // Book total number of events histogram
  st_num_events = new TH1I("st_num_events","ST Number of events; ; Counts",1, 0.5, 1.5);

  //////////////////////// KM hists /////////////////////////////////////
  hnST = new TH1I("hnST",";# of ST digihits;",50,-0.5,49.5);
  hnBCAL = new TH1I("hnBCAL",";# of BCAL digihits;",100,-0.5,99.5);
  hnFCAL = new TH1I("hnFCAL",";# of FCAL digihits;",50,-0.5,49.5);
  hnCDC = new TH1I("hnCDC",";# of CDC digihits;",300,-0.5,299.5);
  hnFDCcathode = new TH1I("hnFDCcathode",";# of FDCcathode digihits;",1000,-0.5,999.5);
  hnFDCanode = new TH1I("hnFDCanode",";# of FDCanode digihits;",1000,-0.5,999.5);
  hnTOF = new TH1I("hnTOF",";# of TOF digihits;",100,-0.5,99.5);
  
  hnST_TDC = new TH1I("hnST_TDC",";# of ST_TDC digihits;",50,-0.5,49.5);
  hnBCAL_TDC = new TH1I("hnBCAL_TDC",";# of BCAL_TDC digihits;",100,-0.5,99.5);
  
  hnTOF_TDC = new TH1I("hnTOF_TDC",";# of TOF_TDC digihits;",50,-0.5,49.5);

  hnEventInfo = new TH1I("hnEventInfo",";# of EventInfo;",100,-0.5,99.5);
  hnROCInfo = new TH1I("hnROCInfo",";# of ROCInfo;",100,-0.5,99.5);
  hnROCId = new TH1I("hnROCId",";counts for each ROC ID;",100,-0.5,99.5);
  // ST vs BCAL
  hADC_ST_BCAL      = new TH2I("hADC_ST_BCAL","time of ST vs BCAL ADC digihits for both up and down stream;BCAL_{ADC};ST_{ADC}",400,0,400,400,0,400);
  hADC_ST_BCAL_up   = new TH2I("hADC_ST_BCAL_up","time of ST vs BCAL ADC digihits for up and stream;BCAL_{ADC};ST_{ADC}",400,0,400,400,0,400);
  hADC_ST_BCAL_down = new TH2I("hADC_ST_BCAL_down","time of ST vs BCAL ADC digihits for both down stream;BCAL_{ADC};ST_{ADC}",400,0,400,400,0,400);

  hTDC_ST_BCAL      = new TH2I("hTDC_ST_BCAL","time of ST vs BCAL TDC digihits for both up and down stream;BCAL_{TDC};ST_{TDC}",400,0,4800,400,0,4800);
  hTDC_ST_BCAL_up   = new TH2I("hTDC_ST_BCAL_up","time of ST vs BCAL TDC digihits for up and stream;BCAL_{TDC};ST_{TDC}",400,0,4800,400,0,4800);
  hTDC_ST_BCAL_down = new TH2I("hTDC_ST_BCAL_down","time of ST vs BCAL TDC digihits for both down stream;BCAL_{TDC};ST_{TDC}",400,0,4800,400,0,4800);

  // ST vs FCAL
  hADC_ST_FCAL      = new TH2I("hADC_ST_FCAL","time of ST vs FCAL ADC digihits for both up and down stream;FCAL_{ADC};ST_{ADC}",400,0,400,400,0,400);

  // ST vs TOF
  hADC_ST_TOF      = new TH2I("hADC_ST_TOF","time of ST vs TOF ADC digihits for both up and down stream;TOF_{ADC};ST_{ADC}",400,0,400,400,0,400);
  hTDC_ST_TOF      = new TH2I("hTDC_ST_TOF","time of ST vs TOF TDC digihits for both up and down stream;TOF_{TDC};ST_{TDC}",400,0,4800,400,0,4800);

  // ST ADC vs ST TDC
  hST_ADC_TDC = new TH2I("hST_ADC_TDC","ST ADC vs TDC;ST_{TDC};ST_{ADC}",400,0,4800,400,0,400);

  gEventInfoAvgTime = new TGraph();
  gEventInfoAvgTime->SetName("gEventInfoAvgTime");

  char gname[200];
  for(Int_t i=0;i<96;i++){
    gtimestamp[i] = new TGraph();
    sprintf(gname,"gtimestamp%2.2d",i);
    gtimestamp[i]->SetName(gname);
  }

  // cd back to main directory
  main->cd();

  outtree = new TTree("timingTree","timingTree");

  // # of hits
  outtree->Branch("nST",&nST,"nST/I");
  outtree->Branch("nBCAL",&nBCAL,"nBCAL/I");
  outtree->Branch("nFCAL",&nFCAL,"nFCAL/I");
  // outtree->Branch("nCDC",&nCDC,"nCDC/I");
  // outtree->Branch("nFDCcathode",&nFDCcathode,"nFDCcathode/I");
  // outtree->Branch("nFDCanode",&nFDCanode,"nFDCanode/I");
  outtree->Branch("nTOF",&nTOF,"nTOF/I");
  outtree->Branch("nTAGH",&nTAGH,"nTAGH/I");
  outtree->Branch("nTAGM",&nTAGM,"nTAGM/I");
  outtree->Branch("nPS",&nPS,"nPS/I");
  outtree->Branch("nPSC",&nPSC,"nPSC/I");

  outtree->Branch("nST_TDC",&nST_TDC,"nST_TDC/I");
  outtree->Branch("nBCAL_TDC",&nBCAL_TDC,"nBCAL_TDC/I");
  outtree->Branch("nTOF_TDC",&nTOF_TDC,"nTOF_TDC/I");
  outtree->Branch("nTAGH_TDC",&nTAGH_TDC,"nTAGH_TDC/I");
  outtree->Branch("nTAGM_TDC",&nTAGM_TDC,"nTAGM_TDC/I");
  outtree->Branch("nPSC_TDC",&nPSC_TDC,"nPSC_TDC/I");

  // ADC times
  outtree->Branch("st_adc_time",st_adc_time,"st_adc_time[nST]/F");
  outtree->Branch("bcal_adc_time",bcal_adc_time,"bcal_adc_time[nBCAL]/F");
  outtree->Branch("fcal_adc_time",fcal_adc_time,"fcal_adc_time[nFCAL]/F");
  outtree->Branch("tof_adc_time",tof_adc_time,"tof_adc_time[nTOF]/F");
  // outtree->Branch("cdc_adc_time",cdc_adc_time,"cdc_adc_time[nCDC]/F");
  // outtree->Branch("fdccathode_adc_time",fdccathode_adc_time,"fdccathode_adc_time[nFDCcathode]/F");
  // outtree->Branch("fdcanode_adc_time",fdcanode_adc_time,"fdcanode_adc_time[nFDCanode]/F");
  outtree->Branch("tagh_adc_time",tagh_adc_time,"tagh_adc_time[nTAGH]/F");
  outtree->Branch("tagm_adc_time",tagm_adc_time,"tagm_adc_time[nTAGM]/F");
  outtree->Branch("ps_adc_time",ps_adc_time,"ps_adc_time[nPS]/F");
  outtree->Branch("psc_adc_time",psc_adc_time,"psc_adc_time[nPSC]/F");

  // TDC times
  outtree->Branch("st_tdc_time",st_tdc_time,"st_tdc_time[nST_TDC]/F");
  outtree->Branch("bcal_tdc_time",bcal_tdc_time,"bcal_tdc_time[nBCAL_TDC]/F");
  outtree->Branch("tof_tdc_time",tof_tdc_time,"tof_tdc_time[nTOF_TDC]/F");
  outtree->Branch("tagh_tdc_time",tagh_tdc_time,"tagh_tdc_time[nTAGH_TDC]/F");
  outtree->Branch("tagm_tdc_time",tagm_tdc_time,"tagm_tdc_time[nTAGM_TDC]/F");
  outtree->Branch("psc_tdc_time",psc_tdc_time,"psc_tdc_time[nPSC_TDC]/F");

  // pulse_integral
  outtree->Branch("st_adc_pulse_integral",st_adc_pulse_integral,"st_adc_pulse_integral[nST]/F");
  outtree->Branch("bcal_adc_pulse_integral",bcal_adc_pulse_integral,"bcal_adc_pulse_integral[nBCAL]/F");
  outtree->Branch("tof_adc_pulse_integral",tof_adc_pulse_integral,"tof_adc_pulse_integral[nTOF]/F");
  outtree->Branch("fcal_adc_pulse_integral",fcal_adc_pulse_integral,"fcal_adc_pulse_integral[nFCAL]/F");
  // outtree->Branch("cdc_adc_pulse_integral",cdc_adc_pulse_integral,"cdc_adc_pulse_integral[nCDC]/F");
  // outtree->Branch("fdccathode_adc_pulse_integral",fdccathode_adc_pulse_integral,"fdccathode_adc_pulse_integral[nCDC]/F");
  outtree->Branch("tagh_adc_pulse_integral",tagh_adc_pulse_integral,"tagh_adc_pulse_integral[nTAGH]/F");
  outtree->Branch("tagm_adc_pulse_integral",tagm_adc_pulse_integral,"tagm_adc_pulse_integral[nTAGM]/F");
  outtree->Branch("ps_adc_pulse_integral",ps_adc_pulse_integral,"ps_adc_pulse_integral[nPS]/F");
  outtree->Branch("psc_adc_pulse_integral",psc_adc_pulse_integral,"psc_adc_pulse_integral[nPSC]/F");

  // pedestal
  outtree->Branch("st_adc_pedestal",st_adc_pedestal,"st_adc_pedestal[nST]/F");
  outtree->Branch("bcal_adc_pedestal",bcal_adc_pedestal,"bcal_adc_pedestal[nBCAL]/F");
  outtree->Branch("tof_adc_pedestal",tof_adc_pedestal,"tof_adc_pedestal[nTOF]/F");
  outtree->Branch("fcal_adc_pedestal",fcal_adc_pedestal,"fcal_adc_pedestal[nFCAL]/F");
  // outtree->Branch("cdc_adc_pedestal",cdc_adc_pedestal,"cdc_adc_pedestal[nCDC]/F");
  // outtree->Branch("fdccathode_adc_pedestal",fdccathode_adc_pedestal,"fdccathode_adc_pedestal[nCDC]/F");
  outtree->Branch("tagh_adc_pedestal",tagh_adc_pedestal,"tagh_adc_pedestal[nTAGH]/F");
  outtree->Branch("tagm_adc_pedestal",tagm_adc_pedestal,"tagm_adc_pedestal[nTAGM]/F");
  outtree->Branch("ps_adc_pedestal",ps_adc_pedestal,"ps_adc_pedestal[nPS]/F");
  outtree->Branch("psc_adc_pedestal",psc_adc_pedestal,"psc_adc_pedestal[nPSC]/F");

  // nsamples_integral
  outtree->Branch("st_adc_nsamples_integral",st_adc_nsamples_integral,"st_adc_nsamples_integral[nST]/I");
  outtree->Branch("bcal_adc_nsamples_integral",bcal_adc_nsamples_integral,"bcal_adc_nsamples_integral[nBCAL]/I");
  outtree->Branch("tof_adc_nsamples_integral",tof_adc_nsamples_integral,"tof_adc_nsamples_integral[nTOF]/I");
  outtree->Branch("fcal_adc_nsamples_integral",fcal_adc_nsamples_integral,"fcal_adc_nsamples_integral[nFCAL]/I");
  outtree->Branch("tagh_adc_nsamples_integral",tagh_adc_nsamples_integral,"tagh_adc_nsamples_integral[nTAGH]/I");
  outtree->Branch("tagm_adc_nsamples_integral",tagm_adc_nsamples_integral,"tagm_adc_nsamples_integral[nTAGM]/I");
  outtree->Branch("ps_adc_nsamples_integral",ps_adc_nsamples_integral,"ps_adc_nsamples_integral[nPS]/I");
  outtree->Branch("psc_adc_nsamples_integral",psc_adc_nsamples_integral,"psc_adc_nsamples_integral[nPSC]/I");

  // nsamples_pedestal
  outtree->Branch("st_adc_nsamples_pedestal",st_adc_nsamples_pedestal,"st_adc_nsamples_pedestal[nST]/I");
  outtree->Branch("bcal_adc_nsamples_pedestal",bcal_adc_nsamples_pedestal,"bcal_adc_nsamples_pedestal[nBCAL]/I");
  outtree->Branch("tof_adc_nsamples_pedestal",tof_adc_nsamples_pedestal,"tof_adc_nsamples_pedestal[nTOF]/I");
  outtree->Branch("fcal_adc_nsamples_pedestal",fcal_adc_nsamples_pedestal,"fcal_adc_nsamples_pedestal[nFCAL]/I");
  outtree->Branch("tagh_adc_nsamples_pedestal",tagh_adc_nsamples_pedestal,"tagh_adc_nsamples_pedestal[nTAGH]/I");
  outtree->Branch("tagm_adc_nsamples_pedestal",tagm_adc_nsamples_pedestal,"tagm_adc_nsamples_pedestal[nTAGM]/I");
  outtree->Branch("ps_adc_nsamples_pedestal",ps_adc_nsamples_pedestal,"ps_adc_nsamples_pedestal[nPS]/I");
  outtree->Branch("psc_adc_nsamples_pedestal",psc_adc_nsamples_pedestal,"psc_adc_nsamples_pedestal[nPSC]/I");

  // QF
  outtree->Branch("st_adc_QF",st_adc_QF,"st_adc_QF[nST]/I");
  outtree->Branch("tof_adc_QF",tof_adc_QF,"tof_adc_QF[nTOF]/I");
  outtree->Branch("bcal_adc_QF",bcal_adc_QF,"bcal_adc_QF[nBCAL]/I");
  outtree->Branch("fcal_adc_QF",fcal_adc_QF,"fcal_adc_QF[nFCAL]/I");
  // outtree->Branch("cdc_adc_QF",cdc_adc_QF,"cdc_adc_QF[nCDC]/I");
  // outtree->Branch("fdccathode_adc_QF",fdccathode_adc_QF,"fdccathode_adc_QF[nFDCcathode]/I");
  outtree->Branch("tagh_adc_QF",tagh_adc_QF,"tagh_adc_QF[nTAGH]/I");
  outtree->Branch("tagm_adc_QF",tagm_adc_QF,"tagm_adc_QF[nTAGM]/I");
  outtree->Branch("ps_adc_QF",ps_adc_QF,"ps_adc_QF[nPS]/I");
  outtree->Branch("psc_adc_QF",psc_adc_QF,"psc_adc_QF[nPSC]/I");

  // detector channel info
  // outtree->Branch("fdccathode_channel",fdccathode_channel,"fdccathode_channel[nFDCcathode]/I");
  // outtree->Branch("fdcanode_channel",fdcanode_channel,"fdcanode_channel[nFDCanode]/I");

  // CDC channel info
  // outtree->Branch("cdc_ring",cdc_ring,"cdc_ring[nCDC]/I");
  // outtree->Branch("cdc_straw",cdc_straw,"cdc_straw[nCDC]/I");

  // ST channel info
  outtree->Branch("st_adc_channel",st_adc_channel,"st_adc_channel[nST]/I");
  outtree->Branch("st_tdc_channel",st_tdc_channel,"st_tdc_channel[nST_TDC]/I");

  // BCAL channel info
  outtree->Branch("bcal_adc_layer",bcal_adc_layer,"bcal_adc_layer[nBCAL]/I");
  outtree->Branch("bcal_adc_sector",bcal_adc_sector,"bcal_adc_sector[nBCAL]/I");
  outtree->Branch("bcal_adc_module",bcal_adc_module,"bcal_adc_module[nBCAL]/I");
  outtree->Branch("bcal_adc_end",bcal_adc_end,"bcal_adc_end[nBCAL]/I");
  outtree->Branch("bcal_tdc_layer",bcal_tdc_layer,"bcal_tdc_layer[nBCAL_TDC]/I");
  outtree->Branch("bcal_tdc_sector",bcal_tdc_sector,"bcal_tdc_sector[nBCAL_TDC]/I");
  outtree->Branch("bcal_tdc_module",bcal_tdc_module,"bcal_tdc_module[nBCAL_TDC]/I");
  outtree->Branch("bcal_tdc_end",bcal_tdc_end,"bcal_tdc_end[nBCAL_TDC]/I");

  // TOF channel info
  outtree->Branch("tof_adc_plane",tof_adc_plane,"tof_adc_plane[nTOF]/I");
  outtree->Branch("tof_tdc_plane",tof_tdc_plane,"tof_tdc_plane[nTOF_TDC]/I");
  outtree->Branch("tof_adc_bar",tof_adc_bar,"tof_adc_bar[nTOF]/I");
  outtree->Branch("tof_tdc_bar",tof_tdc_bar,"tof_tdc_bar[nTOF_TDC]/I");
  outtree->Branch("tof_adc_end",tof_adc_end,"tof_adc_end[nTOF]/I");
  outtree->Branch("tof_tdc_end",tof_tdc_end,"tof_tdc_end[nTOF_TDC]/I");

  // raw data for TOF (only when getRawData is set)
  if(getRawData)
    outtree->Branch("waveform_tof_adc",waveform_tof_adc,"waveform_bcal_adc[nTOF][100]/s");

  // FCAL channel info
  outtree->Branch("fcal_channel",fcal_channel,"fcal_channel[nFCAL]/I");
  outtree->Branch("fcal_row",fcal_row,"fcal_row[nFCAL]/I");
  outtree->Branch("fcal_column",fcal_column,"fcal_column[nFCAL]/I");

  // TAGH channel info
  outtree->Branch("tagh_adc_counter_id",tagh_adc_counter_id,"tagh_adc_counter_id[nTAGH]/I");
  outtree->Branch("tagh_tdc_counter_id",tagh_tdc_counter_id,"tagh_tdc_counter_id[nTAGH_TDC]/I");

  // TAGM channel info
  outtree->Branch("tagm_adc_row",tagm_adc_row,"tagm_adc_row[nTAGM]/I");
  outtree->Branch("tagm_adc_column",tagm_adc_column,"tagm_adc_column[nTAGM]/I");
  outtree->Branch("tagm_tdc_row",tagm_tdc_row,"tagm_tdc_row[nTAGM_TDC]/I");
  outtree->Branch("tagm_tdc_column",tagm_tdc_column,"tagm_tdc_column[nTAGM_TDC]/I");

  // PS channel info
  outtree->Branch("ps_adc_arm",ps_adc_arm,"ps_adc_arm[nPS]/I");
  outtree->Branch("ps_adc_column",ps_adc_column,"ps_adc_column[nPS]/I");

  // PSC channel info
  outtree->Branch("psc_adc_id",psc_adc_id,"psc_adc_id[nPSC]/I");
  outtree->Branch("psc_tdc_id",psc_tdc_id,"psc_tdc_id[nPSC_TDC]/I");

  // outtree->Branch("waveform_bcal_adc",waveform_bcal_adc,"waveform_bcal_adc[nBCAL][100]/i");
  // outtree->Branch("waveform_tagh_adc",waveform_tagh_adc,"waveform_tagh_adc[nTAGH][100]/i");

  // f1TDC times
  outtree->Branch("nf1TDC_TD",&nf1TDC_TD,"nf1TDC_TD/I");
  outtree->Branch("f1tdc_TD",f1tdc_TD,"f1tdc_TD[nf1TDC_TD]/l");
  outtree->Branch("nf1TDC_TS",&nf1TDC_TS,"nf1TDC_TS/I");
  outtree->Branch("f1tdc_TS",f1tdc_TS,"f1tdc_TS[nf1TDC_TS]/l");

  // DCODAEventInfo
  outtree->Branch("nDCODAEventInfo",&nDCODAEventInfo,"nDCODAEventInfo/I");
  outtree->Branch("eventInfoAvgTime",&eventInfoAvgTime,"eventInfoAvgTime/l");
  outtree->Branch("eventNum",&eventNum,"eventNum/I");

  // DCODAROCInfo
  outtree->Branch("nDCODAROCInfo",&nDCODAROCInfo,"nDCODAROCInfo/I");
  outtree->Branch("rocid",rocid,"rocid[96]/I");
  outtree->Branch("rocTime",rocTime,"rocTime[96]/l");

  japp->RootUnLock(); //RELEASE ROOT LOCK!!

  cout << "end of init" << endl;
  return NOERROR;
}
//----------------------------------------------------------------------------------
jerror_t JEventProcessor_timing::brun(JEventLoop *eventLoop, int runnumber) {
  // This is called whenever the run number changes
  return NOERROR;
}
//----------------------------------------------------------------------------------
jerror_t JEventProcessor_timing::evnt(JEventLoop *eventLoop, int eventnumber) {

  // Get all data objects first so we minimize the time we hold the ROOT mutex lock
  // cout << "------------- entry = " << eventnumber << " -------------" << endl;

  // Each detector's hits
  vector<const DSCDigiHit*>         dscdigihits;         // ST DigiHits
  vector<const DSCTDCDigiHit*>      dsctdcdigihits;      // ST TDC DigiHits
  vector<const DSCHit*>             dschits;             // ST

  vector<const DBCALDigiHit*>       dbcaldigihits;       // BCAL DigiHits
  vector<const DBCALTDCDigiHit*>    dbcaltdcdigihits;    // BCAL TDC DigiHits

  vector<const DFCALDigiHit*>       dfcaldigihits;       // FCAL DigiHits

  vector<const DCDCDigiHit*>        dcdcdigihits;        // CDC DigiHits

  vector<const DFDCWireDigiHit*>    dfdcanodedigihits;   // FDC cathode DigiHits
  vector<const DFDCCathodeDigiHit*> dfdccathodedigihits; // FDC anode DigiHits

  vector<const DTOFDigiHit*>        dtofdigihits;        // TOF DigiHits
  vector<const DTOFTDCDigiHit*>     dtoftdcdigihits;     // TOF TDC DigiHits

  vector<const DTAGHDigiHit*>       dtaghdigihits;       // TAGH DigiHits
  vector<const DTAGHTDCDigiHit*>    dtaghtdcdigihits;    // TAGHTDC DigiHits

  vector<const DTAGMDigiHit*>       dtagmdigihits;       // TAGM DigiHits
  vector<const DTAGMTDCDigiHit*>    dtagmtdcdigihits;    // TAGMTDC DigiHits

  vector<const DPSCDigiHit*>        dpscdigihits;        // PSC DigiHits
  vector<const DPSCTDCDigiHit*>     dpsctdcdigihits;     // PSC TDC DigiHits

  vector<const DPSDigiHit*>         dpsdigihits;         // PS DigiHits

  vector<const DCODAEventInfo*>     dcodaeventinfo;      // DCODAEventInfo
  vector<const DCODAROCInfo*>       dcodarocinfo;        // DCODAROCInfo

  // f1TDC hits for timing cable
  vector<const DF1TDCHit *>f1tdchit;

  // waveforms for ADCs
  const Df250PulseIntegral *pulseintegral;
  const Df250WindowRawData *windowrawdata;

  // ST
  eventLoop->Get(dscdigihits);
  eventLoop->Get(dsctdcdigihits);
  eventLoop->Get(dschits);
  // BCAL
  eventLoop->Get(dbcaldigihits);
  eventLoop->Get(dbcaltdcdigihits);
  // FCAL
  eventLoop->Get(dfcaldigihits);
  // CDC
  eventLoop->Get(dcdcdigihits);
  // FDC
  eventLoop->Get(dfdccathodedigihits);
  eventLoop->Get(dfdcanodedigihits);
  // TOF
  eventLoop->Get(dtofdigihits);
  eventLoop->Get(dtoftdcdigihits);

  // PSC
  eventLoop->Get(dpscdigihits);
  eventLoop->Get(dpsctdcdigihits);

  if(dpscdigihits.size()>0) cout << "AAAAA" << endl;
  if(dpsctdcdigihits.size()>0) cout << "BBBBB" << endl;

  // PS
  eventLoop->Get(dpsdigihits);
  if(dpsdigihits.size()>0) cout << "CCCCCC" << endl;

  // TAGH
  eventLoop->Get(dtaghdigihits);
  eventLoop->Get(dtaghtdcdigihits);

  // TAGM
  eventLoop->Get(dtagmdigihits);
  eventLoop->Get(dtagmtdcdigihits);

  // f1TDC hits
  eventLoop->Get(f1tdchit);
  
  eventLoop->Get(dcodaeventinfo);
  eventLoop->Get(dcodarocinfo);

  // Lock ROOT mutex so other threads won't interfere 
  japp->RootWriteLock();
  {

  nf1TDC_TD = 0;
  nf1TDC_TS = 0;

  for(unsigned int i=0;i<f1tdchit.size();i++){
    if(f1tdchit[i]->rocid == 51
       &&
       f1tdchit[i]->slot == 17
       &&
       f1tdchit[i]->channel == 8){
      f1tdc_TD[nf1TDC_TD] = (ULong64_t)f1tdchit[i]->time; // *0.0559; // in ns
      // cout << "event = " << eventnumber << ", channel = 9, f1tdc TD time = " << f1tdc_TD[nf1TDC] << endl;
      nf1TDC_TD++;
    }

    if(f1tdchit[i]->rocid == 51
       &&
       f1tdchit[i]->slot == 17
       &&
       f1tdchit[i]->channel == 15){
      f1tdc_TS[nf1TDC_TS] = (ULong64_t)f1tdchit[i]->time; // *0.0559; // in ns
      // f1tdc_time[nf1TDC] = (ULong64_t)f1tdchit[i]->time; // *0.0559; // in ns
      // cout << "event = " << eventnumber << ", channel = 15, f1tdc time = " << f1tdchit[i]->time << endl;
      nf1TDC_TS++;
    }
  }
  // cout << "nf1TDC_TD = " << nf1TDC_TD << endl;
  // cout << "nf1TDC_TS = " << nf1TDC_TS << endl;

  nDCODAEventInfo = dcodaeventinfo.size();
  hnEventInfo->Fill(nDCODAEventInfo);
  assert(dcodaeventinfo.size()<=1);
  for(UInt_t i=0;i<dcodaeventinfo.size();i++){
    // assert(dcodaeventinfo[i]->run_number==0);
    // assert(dcodaeventinfo[i]->run_type==0);
    // assert(dcodaeventinfo[i]->event_type==0);
    // cout << "run num    = " << dcodaeventinfo[i]->run_number << endl;
    // cout << "run type   = " << dcodaeventinfo[i]->run_type << endl;
    // cout << "event num  = " << dcodaeventinfo[i]->event_number << endl;
    // cout << "event type = " << dcodaeventinfo[i]->event_type << endl;
    // cout << "avg time   = " << dcodaeventinfo[i]->avg_timestamp << endl;
    eventInfoAvgTime = dcodaeventinfo[i]->avg_timestamp;
    eventNum = dcodaeventinfo[i]->event_number;
    gEventInfoAvgTime->SetPoint(gEventInfoAvgTime->GetN(),dcodaeventinfo[i]->event_number,eventInfoAvgTime);
  }
  nDCODAROCInfo = dcodarocinfo.size();
  hnROCInfo->Fill(nDCODAROCInfo);
  // cout << "size of dcodarocinfo: " << nDCODAROCInfo << endl;

  // If DCODAEventInfo is not there, assert that DCODAROCInfo is also not
  if(nDCODAEventInfo==0) assert(nDCODAROCInfo==0);

  for(UInt_t i=0;i<dcodarocinfo.size();i++){
    rocid[i] = dcodarocinfo[i]->rocid;
    hnROCId->Fill(rocid[i]);
    rocTime[i] = dcodarocinfo[i]->timestamp;

    // size of misc is always 0
    assert(dcodarocinfo[i]->misc.size()==0);

    // cout << "rocid     = " << dcodarocinfo[i]->rocid << endl;
    // cout << "timestamp = " << dcodarocinfo[i]->timestamp << endl;
    // cout << "Nmisc     = " << dcodarocinfo[i]->misc.size() << endl;
    // if(dcodarocinfo[0]->timestamp!=0 && dcodarocinfo[i]->timestamp!=0)
    // assert(dcodarocinfo[0]->timestamp == dcodarocinfo[i]->timestamp);
    gtimestamp[rocid[i]]->SetPoint(gtimestamp[rocid[i]]->GetN(),dcodaeventinfo[0]->event_number,dcodarocinfo[i]->timestamp);
  }

  nST = dscdigihits.size();
  nBCAL = dbcaldigihits.size();
  nFCAL = dfcaldigihits.size();
  nCDC = dcdcdigihits.size();
  nFDCcathode = dfdccathodedigihits.size();
  nFDCanode = dfdcanodedigihits.size();
  nTOF = dtofdigihits.size();
  nTAGH = dtaghdigihits.size();
  nTAGM = dtagmdigihits.size();
  nPS = dpsdigihits.size();
  nPSC = dpscdigihits.size();

  nST_TDC = dsctdcdigihits.size();
  nBCAL_TDC = dbcaltdcdigihits.size();
  nTOF_TDC = dtoftdcdigihits.size();
  nTAGH_TDC = dtaghtdcdigihits.size();
  nTAGM_TDC = dtagmtdcdigihits.size();
  nPSC_TDC = dpsctdcdigihits.size();

  // make sure we are within limits of arrays
  assert(nST < NMAX);
  assert(nBCAL < NMAX);
  assert(nFCAL < NMAX);
  assert(nCDC < NMAX_DC);
  assert(nFDCcathode < NMAX_DC);
  assert(nFDCanode < NMAX_DC);
  assert(nTAGH < NMAX);
  assert(nTAGM < NMAX);
  assert(nPS < NMAX);
  assert(nPSC < NMAX);

  assert(nST_TDC < NMAX);
  assert(nBCAL_TDC < NMAX);
  assert(nTOF_TDC < NMAX);
  assert(nTAGH_TDC < NMAX);
  assert(nTAGM_TDC < NMAX);
  assert(nPSC_TDC < NMAX);

  //////////////////////// KM hists /////////////////////////////////////
  hnST->Fill(nST);
  hnBCAL->Fill(nBCAL);
  hnFCAL->Fill(nFCAL);
  hnCDC->Fill(nCDC);
  hnFDCcathode->Fill(nFDCcathode);
  hnFDCanode->Fill(nFDCanode);
  hnTOF->Fill(nTOF);
  
  hnST_TDC->Fill(nST_TDC);
  hnBCAL_TDC->Fill(nBCAL_TDC);
  
  hnTOF_TDC->Fill(nTOF_TDC);

  // ST ADC hits
  for(Int_t st = 0;st<nST;st++){
    st_adc_time[st] = dscdigihits[st]->pulse_time;
    st_adc_channel[st] = dscdigihits[st]->sector;
    st_adc_pulse_integral[st] = dscdigihits[st]->pulse_integral;
    st_adc_pedestal[st] = dscdigihits[st]->pedestal;
    st_adc_nsamples_integral[st] = dscdigihits[st]->nsamples_integral;
    st_adc_nsamples_pedestal[st] = dscdigihits[st]->nsamples_pedestal;
    st_adc_QF[st] = dscdigihits[st]->QF;
  }

  // BCAL ADC hits
  for(Int_t bcal = 0;bcal<nBCAL;bcal++){
    bcal_adc_time[bcal]   = dbcaldigihits[bcal]->pulse_time;
    bcal_adc_layer[bcal]  = dbcaldigihits[bcal]->layer;
    bcal_adc_sector[bcal] = dbcaldigihits[bcal]->sector;
    bcal_adc_module[bcal] = dbcaldigihits[bcal]->module;
    bcal_adc_end[bcal]    = dbcaldigihits[bcal]->end;

    bcal_adc_pulse_integral[bcal] = dbcaldigihits[bcal]->pulse_integral;
    bcal_adc_pedestal[bcal] = dbcaldigihits[bcal]->pedestal;
    bcal_adc_nsamples_integral[bcal] = dbcaldigihits[bcal]->nsamples_integral;
    bcal_adc_nsamples_pedestal[bcal] = dbcaldigihits[bcal]->nsamples_pedestal;
    bcal_adc_QF[bcal] = dbcaldigihits[bcal]->QF;

    // Get waveform
    // (this part of the code is copied from programs/Analysis/plugins/DAQTreeBCAL)
    dbcaldigihits[bcal]->GetSingle(pulseintegral);
    pulseintegral->GetSingle(windowrawdata);

    // Get a vector of the samples for this channel
    const vector<uint16_t> &samplesvector = windowrawdata->samples;
    assert(samplesvector.size() == 100);
    for(UInt_t i=0;i<samplesvector.size();i++){
      // waveform_bcal_adc[bcal][i] = samplesvector[i];
      // cout << "i = " << i << "   " << waveform_bcal_adc[bcal][i] << endl;
    }
  }

  // FCAL ADC hits
  vector< const DFCALGeometry* > geomVec;
  eventLoop->Get( geomVec );
  const DFCALGeometry& fcalGeom = *(geomVec[0]);

  for(Int_t fcal = 0;fcal<nFCAL;fcal++){
    fcal_adc_time[fcal] = dfcaldigihits[fcal]->pulse_time;
    fcal_channel[fcal] = fcalGeom.channel( dfcaldigihits[fcal]->row, dfcaldigihits[fcal]->column );
    fcal_row[fcal] = dfcaldigihits[fcal]->row;
    fcal_column[fcal] = dfcaldigihits[fcal]->column;

    fcal_adc_pulse_integral[fcal] = dfcaldigihits[fcal]->pulse_integral;
    fcal_adc_pedestal[fcal] = dfcaldigihits[fcal]->pedestal;
    fcal_adc_nsamples_integral[fcal] = dfcaldigihits[fcal]->nsamples_integral;
    fcal_adc_nsamples_pedestal[fcal] = dfcaldigihits[fcal]->nsamples_pedestal;
    fcal_adc_QF[fcal] = dfcaldigihits[fcal]->QF;
  }

  // TOF hits
  for(Int_t tof = 0;tof<nTOF;tof++){
    // cout << "processing tof hit " << tof << " / " << nTOF << endl;
    tof_adc_time[tof]  = dtofdigihits[tof]->pulse_time;
    tof_adc_plane[tof] = dtofdigihits[tof]->plane;
    tof_adc_bar[tof]   = dtofdigihits[tof]->bar;
    tof_adc_end[tof]   = dtofdigihits[tof]->end;

    tof_adc_pulse_integral[tof] = dtofdigihits[tof]->pulse_integral;
    tof_adc_pedestal[tof] = dtofdigihits[tof]->pedestal;
    tof_adc_nsamples_integral[tof] = dtofdigihits[tof]->nsamples_integral;
    tof_adc_nsamples_pedestal[tof] = dtofdigihits[tof]->nsamples_pedestal;
    tof_adc_QF[tof] = dtofdigihits[tof]->QF;

    if(getRawData){
      // cout << "getting rawdata for TOF" << endl;
      dtofdigihits[tof]->GetSingle(pulseintegral);
      pulseintegral->GetSingle(windowrawdata);
      // cout << "got rawdata for TOF" << endl;
      // cout << "looping over samples" << endl;
      const vector<uint16_t> &samplesvector = windowrawdata->samples;
      assert(samplesvector.size() == 100);
      for(UInt_t sample=0;sample<samplesvector.size();sample++){
	// cout << "sample = " << sample << " value = " << samplesvector[sample] << endl;
	waveform_tof_adc[tof][sample] = (UShort_t)samplesvector[sample];
	cout << "samplesvector = " << samplesvector[sample] << " saved value = " << waveform_tof_adc[tof][sample] << endl;
      }
    }
  }

  // TAGH ADC hits
  for(Int_t tagh = 0;tagh<nTAGH;tagh++){

    // From DTAGHHit_factory
    tagh_adc_time[tagh] = dtaghdigihits[tagh]->pulse_time;
    tagh_adc_pulse_integral[tagh] = dtaghdigihits[tagh]->pulse_integral;
    tagh_adc_pedestal[tagh] = dtaghdigihits[tagh]->pedestal;
    tagh_adc_nsamples_integral[tagh] = dtaghdigihits[tagh]->nsamples_integral;
    tagh_adc_nsamples_pedestal[tagh] = dtaghdigihits[tagh]->nsamples_pedestal;
    tagh_adc_QF[tagh] = dtaghdigihits[tagh]->QF;

    tagh_adc_counter_id[tagh] = dtaghdigihits[tagh]->counter_id;

    // Get waveform
    // (this part of the code is copied from programs/Analysis/plugins/DAQTreeBCAL)
    dtaghdigihits[tagh]->GetSingle(pulseintegral);
    pulseintegral->GetSingle(windowrawdata);

    // Get a vector of the samples for this channel
    const vector<uint16_t> &samplesvector = windowrawdata->samples;
    assert(samplesvector.size() == 100);
    // cout << "------------- entry = " << eventnumber << " -------------" << endl;
    for(UInt_t i=0;i<samplesvector.size();i++){
      // waveform_tagh_adc[tagh][i] = samplesvector[i];
      // cout << "i = " << i << "   " << waveform_tagh_adc[tagh][i] << endl;
    }
  }

  // TAGM ADC hits
  for(Int_t tagm = 0;tagm<nTAGM;tagm++){
    // From DTAGMHit_factory
    tagm_adc_time[tagm] = dtagmdigihits[tagm]->pulse_time;
    tagm_adc_pulse_integral[tagm] = dtagmdigihits[tagm]->pulse_integral;
    tagm_adc_pedestal[tagm] = dtagmdigihits[tagm]->pedestal;
    tagm_adc_nsamples_integral[tagm] = dtagmdigihits[tagm]->nsamples_integral;
    tagm_adc_nsamples_pedestal[tagm] = dtagmdigihits[tagm]->nsamples_pedestal;
    tagm_adc_QF[tagm] = dtagmdigihits[tagm]->QF;

    tagm_adc_row[tagm] = dtagmdigihits[tagm]->row;
    tagm_adc_column[tagm] = dtagmdigihits[tagm]->column;
  }

  // PS ADC hits
  for(Int_t ps = 0;ps<nPS;ps++){
    // From DPSHit_factory
    ps_adc_time[ps] = dpsdigihits[ps]->pulse_time;
    ps_adc_pulse_integral[ps] = dpsdigihits[ps]->pulse_integral;
    ps_adc_pedestal[ps] = dpsdigihits[ps]->pedestal;
    ps_adc_nsamples_integral[ps] = dpsdigihits[ps]->nsamples_integral;
    ps_adc_nsamples_pedestal[ps] = dpsdigihits[ps]->nsamples_pedestal;
    ps_adc_QF[ps] = dpsdigihits[ps]->QF;

    ps_adc_arm[ps] = dpsdigihits[ps]->arm;
    ps_adc_column[ps] = dpsdigihits[ps]->column;
  }

  // PSC ADC hits
  for(Int_t psc = 0;psc<nPSC;psc++){
    // From DPSCHit_factory
    psc_adc_time[psc] = dpscdigihits[psc]->pulse_time;
    psc_adc_pulse_integral[psc] = dpscdigihits[psc]->pulse_integral;
    psc_adc_pedestal[psc] = dpscdigihits[psc]->pedestal;
    psc_adc_nsamples_integral[psc] = dpscdigihits[psc]->nsamples_integral;
    psc_adc_nsamples_pedestal[psc] = dpscdigihits[psc]->nsamples_pedestal;
    psc_adc_QF[psc] = dpscdigihits[psc]->QF;

    psc_adc_id[psc] = dpscdigihits[psc]->id;
  }

  // CDC ADC hits
  for(Int_t cdc = 0;cdc<nCDC;cdc++){
    // From DCDCHit_factory
    cdc_adc_time[cdc] = dcdcdigihits[cdc]->pulse_time;
    // cdc_adc_channel[cdc] = dscdigihits[cdc]->sector;

    cdc_adc_pulse_integral[cdc] = dcdcdigihits[cdc]->pulse_integral;
    cdc_adc_pedestal[cdc] = dcdcdigihits[cdc]->pedestal;
    cdc_adc_QF[cdc] = dcdcdigihits[cdc]->QF;

    cdc_ring[cdc] = dcdcdigihits[cdc]->ring;
    cdc_straw[cdc] = dcdcdigihits[cdc]->straw;
  }

  // FDC cathode ADC hits
  for(Int_t fdc_cathode = 0;fdc_cathode<nFDCcathode;fdc_cathode++){
    // From DFDCHit_factory
    fdccathode_adc_time[fdc_cathode] = dfdccathodedigihits[fdc_cathode]->pulse_time;
    // fdc_adc_channel[fdc_cathode] = dfdccathodedigihits[fdc_cathode]->sector;
    fdccathode_adc_pulse_integral[fdc_cathode] = dfdccathodedigihits[fdc_cathode]->pulse_integral;
    fdccathode_adc_pedestal[fdc_cathode] = dfdccathodedigihits[fdc_cathode]->pedestal;
    fdccathode_adc_QF[fdc_cathode] = dfdccathodedigihits[fdc_cathode]->QF;
  }

  // FDC anode ADC hits
  for(Int_t fdc_anode = 0;fdc_anode<nFDCanode;fdc_anode++){
    // Taken from FDC_online code:
    // TDC_v3BIN_SIZE = 0.1118 // ns/LSB
    // so TDC_v3BIN_SIZE/8 = 140 [ps],
    // roughly twice of 57 [ps] for TDC_RES?
    fdcanode_adc_time[fdc_anode] = dfdcanodedigihits[fdc_anode]->time;
    // fdc_adc_channel[fdc_anode] = dfdcanodedigihits[fdc_anode]->sector;
  }

  // ST TDC hits
  for(Int_t st_TDC = 0;st_TDC<nST_TDC;st_TDC++){
    st_tdc_time[st_TDC] = dsctdcdigihits[st_TDC]->time;
    st_tdc_channel[st_TDC] = dsctdcdigihits[st_TDC]->sector;
  }

  // BCAL TDC hits
  for(Int_t bcal_TDC = 0;bcal_TDC<nBCAL_TDC;bcal_TDC++){
    bcal_tdc_time[bcal_TDC]   = dbcaltdcdigihits[bcal_TDC]->time;
    bcal_tdc_layer[bcal_TDC]  = dbcaltdcdigihits[bcal_TDC]->layer;
    if(bcal_tdc_layer[bcal_TDC]==0){
      cout << "bcal_tdc_layer[bcal_TDC]==0" << endl;
      abort();
    }
    bcal_tdc_sector[bcal_TDC] = dbcaltdcdigihits[bcal_TDC]->sector;
    bcal_tdc_module[bcal_TDC] = dbcaltdcdigihits[bcal_TDC]->module;
    bcal_tdc_end[bcal_TDC]    = dbcaltdcdigihits[bcal_TDC]->end;
  }

  // TOF TDC hits
  for(Int_t tof_TDC = 0;tof_TDC<nTOF_TDC;tof_TDC++){
    tof_tdc_time[tof_TDC]  = dtoftdcdigihits[tof_TDC]->time;
    tof_tdc_plane[tof_TDC] = dtoftdcdigihits[tof_TDC]->plane;
    tof_tdc_bar[tof_TDC]   = dtoftdcdigihits[tof_TDC]->bar;
    tof_tdc_end[tof_TDC]   = dtoftdcdigihits[tof_TDC]->end;
  }

  // TAGH TDC hits
  for(Int_t tagh_TDC = 0;tagh_TDC<nTAGH_TDC;tagh_TDC++){
    // From DTAGHTDCHit_factory
    tagh_tdc_time[tagh_TDC] = dtaghtdcdigihits[tagh_TDC]->time;
    tagh_tdc_counter_id[tagh_TDC] = dtaghtdcdigihits[tagh_TDC]->counter_id;
  }

  // TAGM TDC hits
  for(Int_t tagm_TDC = 0;tagm_TDC<nTAGM_TDC;tagm_TDC++){
    // From DTAGMTDCHit_factory
    tagm_tdc_time[tagm_TDC] = dtagmtdcdigihits[tagm_TDC]->time;
    tagm_tdc_row[tagm_TDC] = dtagmtdcdigihits[tagm_TDC]->row;
    tagm_tdc_column[tagm_TDC] = dtagmtdcdigihits[tagm_TDC]->column;
  }

  // PSC TDC hits
  for(Int_t psc_TDC = 0;psc_TDC<nPSC_TDC;psc_TDC++){
    // From DPSCTDCHit_factory
    psc_tdc_time[psc_TDC] = dpsctdcdigihits[psc_TDC]->time;
    psc_tdc_id[psc_TDC] = dpsctdcdigihits[psc_TDC]->id;
  }

  //////////////////////////////////////
  //                                  //
  //         ADC comparison           //
  //                                  //
  //////////////////////////////////////
  // ST
  for(Int_t st = 0;st<nST;st++){
    // ST vs BCAL
    for(Int_t bcal = 0;bcal<nBCAL;bcal++){
      hADC_ST_BCAL->Fill(bcal_adc_time[bcal],st_adc_time[st]);

      if(dbcaldigihits[bcal]->end == DBCALGeometry::kUpstream){
	hADC_ST_BCAL_up->Fill(bcal_adc_time[bcal],st_adc_time[st]);
      }else if(dbcaldigihits[bcal]->end == DBCALGeometry::kDownstream){
	hADC_ST_BCAL_down->Fill(bcal_adc_time[bcal],st_adc_time[st]);
      }
    }

    // ST vs FCAL
    for(Int_t fcal = 0;fcal<nFCAL;fcal++){
      hADC_ST_FCAL->Fill(fcal_adc_time[fcal],st_adc_time[st]);
    }

    // ST vs TOF
    for(Int_t tof = 0;tof<nTOF;tof++){
      hADC_ST_TOF->Fill(tof_adc_time[tof],st_adc_time[tof]);
    }

  } // end of loop over ST hits
  
  //////////////////////////////////////
  //                                  //
  //         TDC comparison           //
  //                                  //
  //////////////////////////////////////
  for(Int_t st_TDC = 0;st_TDC<nST_TDC;st_TDC++){
    // ST vs BCAL
    for(Int_t bcal = 0;bcal<nBCAL_TDC;bcal++){
      hTDC_ST_BCAL->Fill(bcal_tdc_time[bcal],st_tdc_time[st_TDC]);
      if(dbcaltdcdigihits[bcal]->end == DBCALGeometry::kUpstream){
	hTDC_ST_BCAL_up->Fill(bcal_tdc_time[bcal],st_tdc_time[st_TDC]);
      }else if(dbcaltdcdigihits[bcal]->end == DBCALGeometry::kDownstream){
	hTDC_ST_BCAL_down->Fill(bcal_tdc_time[bcal],st_tdc_time[st_TDC]);
      }
    }

    // ST vs TOF
    for(Int_t tof_TDC = 0;tof_TDC<nTOF_TDC;tof_TDC++){
      hTDC_ST_TOF->Fill(tof_tdc_time[tof_TDC],st_tdc_time[st_TDC]);
    }

  } // end of loop over ST hits

  //////////////////////////////////////
  //                                  //
  //      ADC vs TDC comparison       //
  //                                  //
  //////////////////////////////////////
    // ST ADC vs TDC
  for(Int_t st = 0;st<nST;st++){
    for(Int_t st_TDC = 0;st_TDC<nST_TDC;st_TDC++){
      if(dscdigihits[st]->sector == dsctdcdigihits[st_TDC]->sector){
	hST_ADC_TDC->Fill(st_tdc_time[st_TDC],st_adc_time[st]);
      }
    }
  }

  // Fill tree!!
  outtree->Fill();
  ////////////////////////////////////////////////////////////////////////////////////////

  // Redefine the hit objects
  uint32_t ADC_hits = dscdigihits.size();
  uint32_t TDC_hits = dsctdcdigihits.size();
  //uint32_t Hits     = dschits.size();

  uint32_t BCAL_TDC_hits = dbcaltdcdigihits.size();


  // Fill multiplicity histograms with the total number of hits
  st_adc_multi->Fill(ADC_hits);
  st_tdc_multi->Fill(TDC_hits);
  st_adc_tdc_multi_2d->Fill(TDC_hits, ADC_hits);

  // Histogram the total number of events
  st_num_events->Fill(1);

  // Clear the multiplicity arrays for each event
  memset(adc_multiplicity_counter, 0, sizeof(adc_multiplicity_counter));
  memset(tdc_multiplicity_counter, 0, sizeof(tdc_multiplicity_counter));

  // DSCDigiHits
  for(uint32_t i = 0; i < ADC_hits; i++)
    {
      const DSCDigiHit *adc_dhit = dscdigihits[i];
      // Calculate some useful quantities
      int hit_channel       = adc_dhit->sector - 1;             // channel hit (ranging from 0 to NCHANNELS-1)
      // uint32_t NSB          = adc_dhit->nsamples_pedestal;      // number of samples for the pedestal (NSB=4)
      //uint32_t avg_pedestal = adc_dhit->pedestal/NSB;           // average pedestal (should be around 100 chan)
      uint32_t avg_pedestal = adc_dhit->pedestal;           // average pedestal (should be around 100 chan)
      uint32_t pulse_time   = adc_dhit->pulse_time*ADC_PT_RES;   // converted pulse time to ns
      
      // Fill 1D Histograms
      st_sec_adc_dhit->Fill(adc_dhit->sector);			// sector digihit from fADC250
      st_pi_dhit->Fill(adc_dhit->pulse_integral);		// pulse integral
      st_pt_dhit->Fill(pulse_time);			// pulse time
      st_ped_dhit->Fill(avg_pedestal);          		// pedestal
      st_qf_dhit->Fill(adc_dhit->QF);				// quality factor
      st_nsi_dhit->Fill(adc_dhit->nsamples_integral);		// nos from pulse integral
      st_nsp_dhit->Fill(adc_dhit->nsamples_pedestal);		// nos from pedestal
      // Fill 2D histograms
      st_pi_dhit_2d->Fill(adc_dhit->sector, adc_dhit->pulse_integral);
      st_ped_dhit_2d->Fill(adc_dhit->sector, avg_pedestal);
      st_pt_dhit_2d->Fill(adc_dhit->sector, pulse_time);
      st_pi_pt_dhit_2d->Fill(adc_dhit->pulse_integral, pulse_time);
      st_pi_nsi_dhit_2d->Fill(adc_dhit->nsamples_integral, adc_dhit->pulse_integral);
      // Fill dynamic array of 1D histograms
      st_pi_dhit_chan[hit_channel]->Fill(adc_dhit->pulse_integral);
      st_ped_dhit_chan[hit_channel]->Fill(avg_pedestal);
      
      // Calculate the fADC250 multiplicities for each channel
      adc_multiplicity_counter[hit_channel] += 1;

      // Aquire the TDC DigiHits
      for(uint32_t i = 0; i < TDC_hits; i++)
	{
	  const DSCTDCDigiHit *tdc_dhit = dsctdcdigihits[i];
          float tdc_dhit_time = tdc_dhit->time*TDC_RES;
	  // Fill 2D histogram
	  st_pi_tdc_dhit_2d->Fill(tdc_dhit_time, adc_dhit->pulse_integral);
	  st_pt_tdc_dhit_2d->Fill(tdc_dhit_time, pulse_time);
	  if (adc_dhit->sector == tdc_dhit->sector)
	    st_sec_hit->Fill(adc_dhit->sector);	 // sector hit with adc and tdc
	}
    }


  // DSCTDCDigiHits
  for(uint32_t i = 0; i < TDC_hits; i++)
    {
      const DSCTDCDigiHit *tdc_dhit = dsctdcdigihits[i];
      float tdc_dhit_time = tdc_dhit->time*TDC_RES;

      // Obtain the average of all BCAL TDC DigiHit Times
      double bcal_time_sum = 0.0;
      double bcal_time_avg = 0.0;
      for(uint32_t j = 0; j < BCAL_TDC_hits; j++)
	{
	  const DBCALTDCDigiHit *bcal_tdc_dhit = dbcaltdcdigihits[j];
	  float bcal_tdc_dhit_time = bcal_tdc_dhit->time*TDC_RES;
	  st_bcal_tdc->Fill(bcal_tdc_dhit_time - tdc_dhit_time);
	  st_bcal_tdc_2d->Fill(tdc_dhit->sector, bcal_tdc_dhit_time - tdc_dhit_time);
	  bcal_time_sum += bcal_tdc_dhit_time;
	  // BCAL upstream hits
	  if(bcal_tdc_dhit->end == DBCALGeometry::kUpstream)
	    {
	      float bcal_tdc_dhit_time_us = bcal_tdc_dhit->time*TDC_RES;
	      st_bcal_tdc_us->Fill(bcal_tdc_dhit_time_us - tdc_dhit_time);
	      st_bcal_tdc_us_2d->Fill(tdc_dhit->sector, bcal_tdc_dhit_time_us - tdc_dhit_time);
	    }
	  //BCAL downstream hits
	  if(bcal_tdc_dhit->end == DBCALGeometry::kDownstream)
	    {
	      float bcal_tdc_dhit_time_ds = bcal_tdc_dhit->time*TDC_RES;
	      st_bcal_tdc_ds->Fill(bcal_tdc_dhit_time_ds - tdc_dhit_time);
	      st_bcal_tdc_ds_2d->Fill(tdc_dhit->sector, bcal_tdc_dhit_time_ds - tdc_dhit_time);
	    }
	}
      // Fill ST/BCAL average time diff histograms
      bcal_time_avg = bcal_time_sum/BCAL_TDC_hits;
      st_bcal_tdc_avg->Fill(bcal_time_avg - tdc_dhit_time);
      st_bcal_tdc_avg_2d->Fill(tdc_dhit->sector, bcal_time_avg - tdc_dhit_time);

      // Fill 1D Histogram
      int hit_channel = tdc_dhit->sector - 1;
      st_sec_tdc_dhit->Fill(tdc_dhit->sector);         // sector hit
      st_tdc_dhit->Fill(tdc_dhit_time);               // tdc time (ns)
      // Fill 2D Histograms
      st_tdc_dhit_2d->Fill(tdc_dhit->sector, tdc_dhit_time);
      // Fill dynamic array of 1D histograms
      st_tdc_dhit_chan[hit_channel]->Fill(tdc_dhit_time);
      // Calculate the f1TDC multiplicities for each channel
      tdc_multiplicity_counter[hit_channel] += 1;
    }

  // Fill the 1D and 2D multiplicity histograms for the fADC250 and f1TDC
  for (uint32_t i = 0; i < NCHANNELS; i++)
    {
      int sector = i + 1;
      if (adc_multiplicity_counter[i] > 0)
	{
	  st_adc_multi_2d->Fill(sector, adc_multiplicity_counter[i]);
	  st_adc_multi_chan[sector - 1]->Fill(adc_multiplicity_counter[i]);
	}

      if (tdc_multiplicity_counter[i] > 0) 
	{
	  st_tdc_multi_2d->Fill(sector, tdc_multiplicity_counter[i]);
	  st_tdc_multi_chan[sector - 1]->Fill(tdc_multiplicity_counter[i]);
	}
    }
  //four(uint32_t
  
  // DSCHits
  /*
  for(uint32_t i = 0; i < dschits.size(); i++) 
    {
      const DSCHit *hit = dschits[i];

          // this is skipped for the moment
          // Fill 1D histograms
          //st_sec_hit->Fill(hit->sector);			 // sector hit
	  //st_de_hit->Fill(hit->dE*1000.0); 			 // dE/dx converted to MeV
	  //st_t_hit->Fill(hit->t);				 // time (ns)
	  //st_sig_t_hit->Fill(hit->sigma_t);			 // uncertanty in time (ns)
	  // Fill 2D histograms
	  //st_t_hit_2d->Fill(hit->sector, hit->t);
	  //st_sig_t_hit_2d->Fill(hit->sector, hit->sigma_t);
	  //st_de_hit_2d->Fill(hit->sector, hit->dE*1000.0);
	  //st_de_t_hit_2d->Fill(hit->t, hit->dE*1000.0);
	  // Fill dynamic array of 1D histograms
	  //st_t_hit_chan[hit->sector - 1]->Fill(hit->t);
	  //st_sig_t_hit_chan[hit->sector - 1]->Fill(hit->sigma_t);
    }
  */
  // Find a hit in two seperate paddles
  if ((TDC_hits == 2) && (ADC_hits == 2))
    {
	  int32_t tdiff_int = dsctdcdigihits[0]->time - dsctdcdigihits[1]->time;
	  float_t tdiff = float(tdiff_int)*TDC_RES;
	  st_tdc_tdiff->Fill(tdiff);
    }
  
  // Find all possible time differences for multiple TDC hits
  if ((TDC_hits > 1) && (ADC_hits > 1) )
    {
      for (uint32_t i = 1; i < TDC_hits ; i++){
	for (uint32_t j = i; j < TDC_hits; j++){
	  int32_t tdiff_int = dsctdcdigihits[j]->time - dsctdcdigihits[i-1]->time;
	  float_t tdiff = float(tdiff_int)*TDC_RES;
	  st_tdc_tdiff_all->Fill(tdiff);
	}
      }
    }

  // Lock ROOT mutex so other threads won't interfere 
  }
  japp->RootUnLock();

  return NOERROR;
}
//----------------------------------------------------------------------------------
jerror_t JEventProcessor_timing::erun(void) {
  // This is called whenever the run number changes, before it is
  // changed to give you a chance to clean up before processing
  // events from the next run number.
  return NOERROR;
}
//----------------------------------------------------------------------------------
jerror_t JEventProcessor_timing::fini(void) {
  gEventInfoAvgTime->Write();
  for(Int_t i=0;i<96;i++){
    // gtimestamp[i]->Write();
  }

  // Called before program exit after event processing is finished.
  return NOERROR;
}
//----------------------------------------------------------------------------------
//----------------------------------------------------------------------------------