// $Id$
//
//    File: JEventProcessor_trig_rate_thresh.cc
// Created: Sat Sep 10 11:32:36 EDT 2016
// Creator: davidl (on Linux gluon48.jlab.org 2.6.32-431.20.3.el6.x86_64 x86_64)
//

#include <iostream>
using namespace std;

#include "JEventProcessor_trig_rate_thresh.h"
using namespace jana;

#include <DAQ/Df250PulseData.h>
#include <DAQ/DParsedEvent.h>
#include <DAQ/DCODAEventInfo.h>
#include <DAQ/DBeamCurrent.h>
#include <TRIGGER/DL1Trigger.h>
#include <FCAL/DFCALDigiHit.h>
#include <BCAL/DBCALDigiHit.h>
#include <FCAL/DFCALHit.h>
#include <BCAL/DBCALHit.h>


// Routine used to create our JEventProcessor
#include <JANA/JApplication.h>
#include <JANA/JFactory.h>
extern "C"{
void InitPlugin(JApplication *app){
	InitJANAPlugin(app);
	app->AddProcessor(new JEventProcessor_trig_rate_thresh());
}
} // "C"


//------------------
// JEventProcessor_trig_rate_thresh (Constructor)
//------------------
JEventProcessor_trig_rate_thresh::JEventProcessor_trig_rate_thresh()
{

}

//------------------
// ~JEventProcessor_trig_rate_thresh (Destructor)
//------------------
JEventProcessor_trig_rate_thresh::~JEventProcessor_trig_rate_thresh()
{

}

//------------------
// init
//------------------
jerror_t JEventProcessor_trig_rate_thresh::init(void)
{
	Ngoodtrigs = 0;
	Ntot = 0;
	
	timestamp_min = 0.0;
	timestamp_max = 0.0;

	// This was apparently used in Spring 2017 as well as in Alex S.
	// trigger studies in Fall 2016 runs 22102-22115. The high intensity
	// testing I did just before that (runs 22062-22068) used the trigger
	// that was default at that time. Empirically:
	// 25*fcal_tot_en + bcal_tot_en > 45000
	TRIG_FCAL_SCALE = 25;
	TRIG_THRESH     = 45000;

	gPARMS->SetDefaultParameter("TRIG_FCAL_SCALE", TRIG_FCAL_SCALE, "Scale FCAL pulse integral sum by this before adding to BCAL sum and applying threshold");
	gPARMS->SetDefaultParameter("TRIG_THRESH", TRIG_THRESH, "Threshold to apply to BCAL+FCAL*TRIG_FCAL_SCALE for trigger simulation");

	sum_pulse_int_fcal = new TH1D("sum_pulse_int_fcal", "Sum of FCAL pulse_int", 2000, -3000.0, 40000.0);
	sum_pulse_int_bcal = new TH1D("sum_pulse_int_bcal", "Sum of BCAL pulse_int", 2000, -24000.0, 330000.0);
	sum_pulse_int_bcal_vs_fcal = new TH2D("sum_pulse_int_bcal_vs_fcal", "BCAL vs. FCAL with cut on L1 trig emulation", 300, 0.0, 15000.0, 300, 0.0, 75000.0);
	sum_pulse_int_bcal_vs_fcal_nocut = new TH2D("sum_pulse_int_bcal_vs_fcal_nocut", "BCAL vs. FCAL without cut on L1 trig emulation", 300, 0.0, 15000.0, 300, 0.0, 75000.0);
	sum_pulse_int_bcal_vs_fcal_trig6 = new TH2D("sum_pulse_int_bcal_vs_fcal_trig6", "BCAL vs. FCAL for trig 6 without cut", 300, 0.0, 15000.0, 300, 0.0, 75000.0);

	Efcal_trig_nocut = new TH1D("Efcal_trig_nocut", "" , 450, -0.5, 8.0);
	Ebcal_trig_nocut = new TH1D("Ebcal_trig_nocut", "" , 450, -0.5, 8.0);
	Efcal_trig = new TH1D("Efcal_trig", "" , 450, -0.5, 8.0);
	Ebcal_trig = new TH1D("Ebcal_trig", "" , 450, -0.5, 8.0);

	daq_event_tdiff = new TH1D("daq_event_tdiff", "Time between events with cut on trig threshold", 10000, 0.0, 1.0E1);
	daq_event_tdiff_nocut = new TH1D("daq_event_tdiff_nocut", "Time between events", 10000, 0.0, 1.0E1);
	
	tt1rate = new TTree("t1rate", "Trig 1 scaler rate from sync events");
	tt1rate->Branch("rate", &t1rate, "rate/F");
	tt1rate->Branch("t",    &t1time, "t/F");
	tt1rate->Branch("bit",  &t1bit,  "bit/I");

	ttdiff = new TTree("ttdiff", "Times of 2 adjacent events");
	ttdiff->Branch("t1", &time1, "t1/F");
	ttdiff->Branch("t2", &time2, "t2/F");

#if 0
	daq_stats = new TH1D("daq_stats", "Some stats", 10, 0.0, 10.0);
	daq_stats->GetXaxis()->SetBinLabel(1, "Nevents");
	daq_stats->GetXaxis()->SetBinLabel(2, "T_tot (s)");
	daq_stats->GetXaxis()->SetBinLabel(3, "T_min (s)");
	daq_stats->GetXaxis()->SetBinLabel(4, "T_max (s)");
	daq_stats->GetXaxis()->SetBinLabel(5, "Avg rate (Hz)");
	
	tevent = new TTree("tevent", "Time difference between 2 adjacent events");
	tevent->Branch("Event", "Event", &evt);
#endif

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_trig_rate_thresh::brun(JEventLoop *eventLoop, int32_t runnumber)
{
	fcal_cell_thr  =  65; // From FCAL config FADC250_TRIG_THR = 165
	bcal_cell_thr  =  20; // From BCAL config FADC250_TRIG_THR = 120

	fcal_row_mask_min = 26;
	fcal_row_mask_max = 32;
	fcal_col_mask_min = 26;
	fcal_col_mask_max = 32;

	return NOERROR;
}

//------------------
// IsFCALinner2Rings
//------------------
bool JEventProcessor_trig_rate_thresh::IsFCALinner2Rings(const DDAQAddress *a)
{
//return false;
	// Returns true if this channel belongs to one of
	// the FCAL 2 inner rings which were excluded from
	// the trigger. Returns false otherwise.

	const uint32_t &rocid = a->rocid;
	const uint32_t &slot  = a->slot;
	const uint32_t &chan  = a->channel;

	// exclude inner two rings from sum
	switch(rocid){
		case 12:
			if(slot==3 && chan>=2 && chan<=15) return true;
			if(slot==5 && (chan==0 || chan==4 || chan==8 || chan==12)) return true;
			break;
		case 15:
			if(slot==3 && chan<=7) return true;
			if(slot==3 && chan>=9 && chan<=11) return true;
			if(slot==3 && chan>=13 && chan<=15) return true;
			if(slot==5 && chan>=12 && chan<=15) return true;
			break;
		case 18:
			if(slot==3 && chan<=13) return true;
			if(slot==5 && (chan==3 || chan==7 || chan==11 || chan==15)) return true;
			break;
		case 21:
			if(slot==3 && chan<=2) return true;
			if(slot==3 && chan>=4 && chan<=6) return true;
			if(slot==3 && chan>=8 && chan<=15) return true;
			if(slot==5 && chan<=3) return true;
			break;
	}

	return false;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_trig_rate_thresh::evnt(JEventLoop *loop, uint64_t eventnumber)
{
	
	vector<const DBCALDigiHit*  > locBCALDigiHits;
	vector<const DFCALDigiHit*  > locFCALDigiHits;
	vector<const DBCALHit*      > locBCALHits;
	vector<const DFCALHit*      > locFCALHits;
	vector<const DCODAEventInfo*> locCODAEventInfos;
	vector<const DL1Trigger*    > locL1Triggers;

	loop->Get(locBCALDigiHits);
	loop->Get(locFCALDigiHits);
	loop->Get(locBCALHits);
	loop->Get(locFCALHits);
	loop->Get(locCODAEventInfos  );
	loop->Get(locL1Triggers  );
	const DL1Trigger* locL1Trigger = locL1Triggers.empty() ? NULL : locL1Triggers[0];	
	if(!locL1Trigger  ) return NOERROR;
	const DCODAEventInfo* codaeventinfo = locCODAEventInfos.empty() ? NULL : locCODAEventInfos[0];	
	if(!codaeventinfo) return NOERROR;
	
	double timestamp = ((double)codaeventinfo->avg_timestamp)/250.0E6;	

	bool trig1 = (locL1Trigger->trig_mask >> 0) & 0x01;
	bool trig6 = (locL1Trigger->trig_mask >> 5) & 0x01;

	// Get total FCAL energy
	int fcal_tot_en = 0;
	for( auto fcal_hit : locFCALDigiHits){
		int row = fcal_hit->row;
		int col = fcal_hit->column;
		if( (row >= fcal_row_mask_min) && (row <= fcal_row_mask_max) ){
			if( (col >= fcal_col_mask_min) && (col <= fcal_col_mask_max) ) continue;
		}

		if( ((int32_t)fcal_hit->pulse_peak-100) <= fcal_cell_thr) continue;

		uint32_t adc_time = (fcal_hit->pulse_time >> 6) & 0x1FF; // consider only course time
		if((adc_time < 15) || (adc_time > 50)) continue; // changed from 20 and 70 based on run 30284  2/5/2017 DL

		Int_t pulse_int = fcal_hit->pulse_integral - fcal_hit->nsamples_integral*100;
		if(pulse_int < 0) continue;
		fcal_tot_en += pulse_int;
	}

	// Get total BCAL energy
	int bcal_tot_en = 0;
	for(auto bcal_hit : locBCALDigiHits){
		if( ((int32_t)bcal_hit->pulse_peak-100) <= bcal_cell_thr) continue;
		Int_t pulse_int = bcal_hit->pulse_integral - bcal_hit->nsamples_integral*100;
		if(pulse_int < 0) continue;
		bcal_tot_en += pulse_int;
	}
	
	// Trig bit 1 with stronger cut (see comment in init method above)
	bool trig1b = TRIG_FCAL_SCALE*fcal_tot_en + bcal_tot_en > TRIG_THRESH;
	
	// Calibrated FCAL hits
	double sum_calib_fcal = 0.0;
	for(auto h : locFCALHits){
		vector<const Df250PulseData*> mypds;
		h->Get(mypds);
		if(mypds.empty()) continue;
		if(IsFCALinner2Rings(mypds[0])) continue;
		sum_calib_fcal += h->E; // convert to GeV
	}

	// Calibrated BCAL hits
	double sum_calib_bcal = 0.0;
	for(auto h : locBCALHits){
		sum_calib_bcal += h->E;
	}

	//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
	japp->RootFillLock(this);
	
	// Event rate
	if(trig1b){
		ts_history.insert(timestamp);
		if(ts_history.size() > 1000){
			double t1 = *(ts_history.begin());
			double t2 = *(ts_history.erase(ts_history.begin()));
			daq_event_tdiff->Fill( (t2-t1)*1.0E3 ); // time difference in ms
		}
	}
	
	t1bit = 1;
	for(auto r : locL1Trigger->gtp_rate){
		if(r!=0) {
			t1rate = r;
			t1time = timestamp;
			tt1rate->Fill();
		}
		t1bit++;
	}

	if(trig1){
		ts_history_nocut.insert(timestamp);
		if(ts_history_nocut.size() > 1000){
			double t1 = *(ts_history_nocut.begin());
			double t2 = *(ts_history_nocut.erase(ts_history_nocut.begin()));
			daq_event_tdiff_nocut->Fill( (t2-t1)*1.0E3 ); // time difference in ms
			
			time1 = t1;
			time2 = t2;
			ttdiff->Fill();
		}
	}

	// Trigger
	sum_pulse_int_fcal->Fill(fcal_tot_en);
	sum_pulse_int_bcal->Fill(bcal_tot_en);
	if(trig1 ) sum_pulse_int_bcal_vs_fcal_nocut->Fill(fcal_tot_en, bcal_tot_en);
	if(trig6 ) sum_pulse_int_bcal_vs_fcal_trig6->Fill(fcal_tot_en, bcal_tot_en);
	if(trig1b) sum_pulse_int_bcal_vs_fcal->Fill(fcal_tot_en, bcal_tot_en);

	Efcal_trig_nocut->Fill(sum_calib_fcal);
	Ebcal_trig_nocut->Fill(sum_calib_bcal);
	if(trig1b) Efcal_trig->Fill(sum_calib_fcal);
	if(trig1b) Ebcal_trig->Fill(sum_calib_bcal);

	japp->RootFillUnLock(this);
	//- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -


#if 0
	vector<const Df250PulsePedestal*> pps;
	vector<const Df250PulseIntegral*> pis;
	vector<const DFCALHit*> fcalhits;
	vector<const DBCALHit*> bcalhits;
	loop->Get(pps);
	loop->Get(pis);
	loop->Get(fcalhits);
	loop->Get(bcalhits);

	const DL1Trigger *l1trig=NULL;
	try{
		loop->GetSingle(l1trig);
	}catch(...){ return NOERROR; }
	if(!l1trig) return NOERROR;
//	if( (l1trig->trig_mask&0x01) != 0x01 ) return NOERROR;
	
	const DCODAEventInfo *codaeventinfo = NULL;
	try{
		loop->GetSingle(codaeventinfo);
	}catch(...){ return NOERROR; }


	// Pulse integral
	double sum_int_fcal = 0;
	double sum_int_bcal = 0;
	double sum_peak_fcal = 0;
	double sum_peak_bcal = 0;
	for(auto pi : pis){
	
		const uint32_t &rocid = pi->rocid;
		const Df250PulseTime *pt = NULL;
		const Df250PulsePedestal *pp = NULL;
		pi->GetSingle(pt);
		pi->GetSingle(pp);
		if(!pt) continue;
		if(!pp) continue;
		
		// FCAL is rocids 11-22
		if(rocid>=11 && rocid<=22){
			
			if(IsFCALinner2Rings(pp)) continue;
			
			// FCAL f250PulseTime hits in time with trigger 
			// are between 3000-3600 (a 37.5ns range)
//			if(pt->time<3000 || pt->time>3600) continue;
			
			sum_peak_fcal += (double)pp->pulse_peak - (double)pp->pedestal;

			double ped = (double)pi->pedestal*(double)pi->nsamples_integral/(double)pi->nsamples_pedestal;
			sum_int_fcal += (double)pi->integral - ped;
		}
		
		// BCAL is rocids 31-41
		if(rocid>=31 && rocid<=41){

			// BCAL f250PulseTime hits in time with trigger
			// are between 2200-3000 (a 50ns range)
//			if(pt->time<2200 || pt->time>3000) continue;
			
			sum_peak_bcal += (double)pp->pulse_peak - (double)pp->pedestal;

			double ped = (double)pi->pedestal*(double)pi->nsamples_integral/(double)pi->nsamples_pedestal;
			sum_int_bcal += (double)pi->integral - ped;
		}		
	}

	// Calibrated FCAL hits
	double sum_calib_fcal = 0.0;
	for(auto h : fcalhits){
		vector<const Df250PulseIntegral*> mypis;
		h->Get(mypis);
		if(pis.empty()) continue;
		if(IsFCALinner2Rings(pis[0])) continue;
		sum_calib_fcal += h->E; // convert to GeV
	}

	// Calibrated BCAL hits
	double sum_calib_bcal = 0.0;
	for(auto h : bcalhits){
		sum_calib_bcal += h->E;
	}
	
	// Simulate trigger cut using pulse integrals
	double trig_amp = sum_int_fcal + TRIG_BCAL_SCALE*sum_int_bcal;
	bool goodtrig = trig_amp >= TRIG_THRESH;

	if(goodtrig) Ngoodtrigs++;
	Ntot++;
	
	// Access DParsedEvent so we can get EVIO buffer info
	void *vptr = loop->GetJEvent().GetRef();
	DParsedEvent *parsedevent = static_cast<DParsedEvent*>(vptr);
	
	
	uint64_t trigbits_lo = l1trig->trig_mask;
	uint64_t trigbits_hi = l1trig->fp_trig_mask;
	uint64_t trigbits = (trigbits_lo) | (trigbits_hi << 32);// 0-31=gtp  32-63=fp
	tsinfo_t tsinfo={l1trig->event_type, trigbits, parsedevent->buff_len, parsedevent->istreamorder};

	japp->RootFillLock(this);



	// Time difference between events
	// since events may be out of order due to L3, we
	// keep track of up to 400 timestamps and only make 
	// entries once we have accumulated that many.
	// (probably better to look for adjacent event numbers
	// but that will take a littel refactoring.)
	const uint64_t &timestamp = codaeventinfo->avg_timestamp;
	ts_history_nocut[timestamp] = tsinfo;
	if(goodtrig) ts_history[timestamp] = tsinfo;
	double tdiff_ms_nocut = 0.0;
	double tdiff_ms = 0.0;
	int event1_type = 0;
	int event2_type = 0;
	uint64_t trig1 = 0;
	uint64_t trig2 = 0;
	uint32_t buff_len1 = 0;
	uint32_t buff_len2 = 0;
	uint64_t block1 = 0;
	uint64_t block2 = 0;
	if(ts_history_nocut.size()>4000){
		auto it1 = ts_history_nocut.begin();
		auto it2 = it1;
		it2++;
		uint64_t t1 = it1->first;
		uint64_t t2 = it2->first;
		event1_type = it1->second.event_type;
		event2_type = it2->second.event_type;
		trig1 = it1->second.trigbits;
		trig2 = it2->second.trigbits;
		buff_len1 = it1->second.buff_len;
		buff_len2 = it2->second.buff_len;
		block1 = it1->second.block;
		block2 = it2->second.block;
		ts_history_nocut.erase(it1, it2);
		double tdiff_ns = (double)(t2 - t1)*4.0;
		tdiff_ms_nocut = tdiff_ns/1.0E6;
	}
	if(ts_history.size()>4000){
		auto it1 = ts_history.begin();
		auto it2 = it1;
		it2++;
		uint64_t t1 = it1->first;
		uint64_t t2 = it2->first;
		ts_history.erase(it1, it2);
		double tdiff_ns = (double)(t2 - t1)*4.0;
		tdiff_ms = tdiff_ns/1.0E6;
	}
	
	// Keep track of min and max timestamps so we can calculate total average
	// event rate for file.
	if( (timestamp_min==0) || (timestamp < timestamp_min) ) timestamp_min = timestamp;
	if( (timestamp_max==0) || (timestamp > timestamp_max) ) timestamp_max = timestamp;
	
	sum_pulse_peak_fcal->Fill(sum_peak_fcal);
	sum_pulse_peak_bcal->Fill(sum_peak_bcal);
	if(goodtrig)sum_pulse_peak_bcal_vs_fcal->Fill(sum_peak_fcal, sum_peak_bcal);
	sum_pulse_peak_bcal_vs_fcal_nocut->Fill(sum_peak_fcal, sum_peak_bcal);

	sum_pulse_int_fcal->Fill(sum_int_fcal);
	sum_pulse_int_bcal->Fill(sum_int_bcal);
	if(goodtrig)sum_pulse_int_bcal_vs_fcal->Fill(sum_int_fcal, sum_int_bcal);
	sum_pulse_int_bcal_vs_fcal_nocut->Fill(sum_int_fcal, sum_int_bcal);
	
	Efcal_trig_nocut->Fill(sum_calib_fcal);
	Ebcal_trig_nocut->Fill(sum_calib_bcal);

	if(goodtrig)Efcal_trig->Fill(sum_calib_fcal);
	if(goodtrig)Ebcal_trig->Fill(sum_calib_bcal);

	if(tdiff_ms_nocut!=0.0) daq_event_tdiff_nocut->Fill(tdiff_ms_nocut);
	if(tdiff_ms!=0.0) daq_event_tdiff->Fill(tdiff_ms);
	
	if(tdiff_ms_nocut!=0.0) {
		evt.tdiff = tdiff_ms_nocut;
		evt.event1_type = event1_type;
		evt.event2_type = event2_type;
		evt.trig1 = trig1;
		evt.trig2 = trig2;
		evt.buff_len1 = buff_len1;
		evt.buff_len2 = buff_len2;
		evt.is_same_block = (block1 == block2);
		tevent->Fill();
	}
	
	japp->RootFillUnLock(this);

#endif

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_trig_rate_thresh::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;
}

//------------------
// fini
//------------------
jerror_t JEventProcessor_trig_rate_thresh::fini(void)
{
#if 0
	cout << "Ngoodtrigs: " << Ngoodtrigs << endl;
	cout << "      Ntot: " << Ntot << endl;
	
	double tmin_sec = (double)timestamp_min*4.0/1.0E9;
	double tmax_sec = (double)timestamp_max*4.0/1.0E9;
	double t_sec = tmax_sec - tmin_sec;
	daq_stats->SetBinContent(1, (double)Ntot);
	daq_stats->SetBinContent(2, t_sec);
	daq_stats->SetBinContent(3, tmin_sec);
	daq_stats->SetBinContent(4, tmax_sec);
	daq_stats->SetBinContent(5, (double)Ntot/t_sec);
#endif	
	
	return NOERROR;
}