// $Id$
//
//    File: JEventProcessor_fmwpcana.cc
// Created: Thu Jan 18 15:25:47 EST 2018
// Creator: hdmuonops (on Linux gluon118.jlab.org 2.6.32-642.3.1.el6.x86_64 x86_64)
//

#include "JEventProcessor_fmwpcana.h"
using namespace jana;

#include <TF1.h>
#include <TMath.h>

#include <iostream>
using namespace std;


// 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_fmwpcana());
}
} // "C"

#include <DAQ/Df125WindowRawData.h>
#include <DAQ/Df250WindowRawData.h>


//------------------
// JEventProcessor_fmwpcana (Constructor)
//------------------
JEventProcessor_fmwpcana::JEventProcessor_fmwpcana()
{

}


//-------------------
// Make125Histos
//-------------------
void JEventProcessor_fmwpcana::Make125Histos(int Nbins)
{

	cout << "Making f125 histograms with " << Nbins << " bins ...." << endl;

	for(int i=0; i<144; i++){
		char hname[256];
		char title[256];
		sprintf(hname, "fmwpc_wire%03d", i);
		sprintf(title, "FMWPC wire %d;sample (8ns);amplitude (~0.122mV)", i);
		auto h = new TH1I(hname, title, Nbins, 0.0, (float)Nbins);
		h->GetYaxis()->SetRangeUser(0.0, 8000.0);
		hfmwpc.push_back( h );
	}
}

//-------------------
// Make250Histos
//-------------------
void JEventProcessor_fmwpcana::Make250Histos(int Nbins)
{

	cout << "Making f250 histograms with " << Nbins << " bins ...." << endl;

	for(int i=0; i<4; i++){
		char hname[256];
		string title;
		sprintf(hname, "paddles_chan%03d", i);
		switch(i){
			case 0:  title = "PMT downstream, bottom";  break;
			case 1:  title = "PMT upstream, bottom";  break;
			case 2:  title = "PMT downstream, top";  break;
			case 3:  title = "PMT upstream, top";  break;
		}
		auto h = new TH1I(hname, title.c_str(), Nbins, 0.0, (float)Nbins);
		h->SetXTitle("sample (4ns)");
		h->SetYTitle("amplitude");
		h->GetYaxis()->SetRangeUser(0.0, 2500.0);
		hpaddles.push_back( h );
		
		char fname[256];
		sprintf(fname, "fpaddles_chan%03d", i);
		auto f = new TF1(fname, "[0]*TMath::Landau(x,[1],[2], {kTRUE})+[3]");
		f->SetParName(0, "amplitude");
		f->SetParName(1, "time");
		f->SetParName(2, "width");
		f->SetParName(3, "pedestal");
		fpaddles.push_back( f );
	}
}

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

}

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

	htpaddles.push_back(new TH1D("tpaddle0", "Fit time paddle downstream, bottom", 800, 0.0, 800.0));
	htpaddles.push_back(new TH1D("tpaddle1", "Fit time paddle upstream, bottom", 800, 0.0, 800.0));
	htpaddles.push_back(new TH1D("tpaddle2", "Fit time paddle downstream, top", 800, 0.0, 800.0));
	htpaddles.push_back(new TH1D("tpaddle3", "Fit time paddle upstream, top", 800, 0.0, 800.0));

	tdiff_upstream = new TH1D("tdiff_upstream", "#Deltat Upstream paddle", 100, -50.0, 50.0);
	tdiff_downstream = new TH1D("tdiff_downstream", "#Deltat Downstream paddle", 100, -50.0, 50.0);
	tdiff_up_vs_down = new TH2D("tdiff_up_vs_down", "#Deltat Upstream vs. Downstream paddle", 100, -10.0, 10.0, 100, -10.0, 10.0);
	tup_minus_tdown = new TH1D("tup_minus_tdown", "#Deltat Upstream minus Downstream paddle", 200, -20.0, 20.0);
	tdiff_means = new TH1D("tdiff_means", "Difference of mean times Upstream-Downstream", 200, -20.0, 20.0);
	tmean_diffs = new TH1D("tmean_diffs", "Mean of time differences", 200, -20.0, 20.0);

	amp_vs_wire = new TH2D("amp_vs_wire", "Peak amp. vs. wire", 144, 0.0, 144.0, 410, 0.0, 4100.0);
	paddletime_vs_wire = new TH2D("paddletime_vs_wire", "Paddle time vs. wire number", 144, 0.0, 144.0, 200, -20.0, 20.0);

	double ns_per_bin = 8.0;
	double dft_min =  0.2; // min frequency in Mhz
	double dft_max = 30.0; // max frequency in Mhz
	int Ndft_bins = (1.0/dft_min - 1.0/dft_max)*1.0E3/ns_per_bin;
	dft_by_wire = new TH2D("dft_by_wire", "Discrete Fourier Transform of noise;Frequncy (MHz);Wire number", Ndft_bins, dft_min, dft_max, 144, 0.0, 144.0);

	// This is used to cut on peak value as a function of wire number.
	// it was derived by eye looking at cosmic run 85.
	//fthresh = new TF1("fun", "[0]*(1.0-[2]*pow(abs(x-[1])-[1],2.0))", 0.0, 144.0);
	fthresh = new TF1("fun", "[0]*(1.0-[2]*pow(x-[1],2.0))", 0.0, 144.0);
	fthresh->SetParameter(0, 1550.0);
	fthresh->SetParameter(1, 73.0);
	fthresh->SetParameter(2, 1050.0/1550.0/pow(75.0,2.0));

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t JEventProcessor_fmwpcana::brun(JEventLoop *eventLoop, int32_t runnumber)
{
	// This is called whenever the run number changes
	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t JEventProcessor_fmwpcana::evnt(JEventLoop *loop, uint64_t eventnumber)
{
	vector<const Df125WindowRawData*> wrds125;
	vector<const Df250WindowRawData*> wrds250;
	loop->Get(wrds125);
	loop->Get(wrds250);


	// Discrete Fourier transform is CPU intensive. Do this
	// outside of mutex lock.
	map<int, vector<double> > dfts;
	for(auto wrd : wrds125){
	
		uint32_t iwire = (71-wrd->channel) + (wrd->slot==3 ? 0:72);
		int Nsamples = wrd->samples.size();
		int Ndft_bins     = dft_by_wire->GetNbinsX();
		double ns_per_bin = 8;
		double domain   = ns_per_bin*(double)Nsamples;
		for(int ibin=1; ibin<=Ndft_bins; ibin++){

			// n.b. the factor of 1.0E3 below is to convert from MHz to GHz
			double f = dft_by_wire->GetXaxis()->GetBinCenter(ibin);
			double k = domain/(1.0E3/f); // k is number of oscillations in whole x-range for this freq.

			double dphi = TMath::TwoPi()*k/(double)Nsamples;

			double sumc = 0.0;
			double sums = 0.0;
			double theta =0.0;
			for( double v : wrd->samples ){
				theta += dphi;
				sumc += v*cos(theta);
				sums += v*sin(theta);
			}

			double w = sqrt(sumc*sumc + sums*sums);
			dfts[iwire].push_back( w );
		}
	}

	japp->RootFillLock(this);

	// ------------------ f125
	if( hfmwpc.empty() && !wrds125.empty() ){
		Make125Histos( wrds125[0]->samples.size() );	
	}
	
	vector<uint32_t> fmwpchits;
	for(auto wrd : wrds125){
		
		uint32_t iwire = (71-wrd->channel) + (wrd->slot==3 ? 0:72);
		if( iwire >= hfmwpc.size() ) continue;
		auto h = hfmwpc[iwire];
		auto Nbins = h->GetNbinsX();
		
		double max = 0.0;
		for(uint32_t i=0; i<wrd->samples.size(); i++){
			if( (int)i<Nbins ) h->SetBinContent(i+1, wrd->samples[i]);
			if( wrd->samples[i]>max ) max = wrd->samples[i];
		}
		if(max>fthresh->Eval(iwire)) fmwpchits.push_back(iwire);
		if(max>1.0) amp_vs_wire->Fill(iwire, max);

	}
	
	// Fill DFT histo
	for(auto p : dfts){
		int iwire = p.first;
		vector<double> &vdft = p.second;
		int ibin=1;
		for( auto w : vdft ){
			double f = dft_by_wire->GetXaxis()->GetBinCenter(ibin++);
			dft_by_wire->Fill(f, iwire, w);
		}
	}

	// ------------------ f250
	if( hpaddles.empty() && !wrds250.empty() ){
		Make250Histos( wrds250[0]->samples.size() );
	}
	
	vector<double> thit = {-1.0, -1.0, -1.0, -1.0};
	for(auto wrd : wrds250){

		uint32_t ichan = wrd->channel;
		if( ichan >= hpaddles.size() ) continue;
		auto h = hpaddles[ichan];
		auto Nbins = h->GetNbinsX();
		
		double max = 0;
		double xmax = 0;
		for(uint32_t i=0; i<wrd->samples.size(); i++){
			if( (int)i<Nbins ) {
				auto val = wrd->samples[i];
				h->SetBinContent(i+1, val);

				if((double)val > max) {
					max = (double)val;
					xmax = (double)i;
				}
			}
		}
		
		if( max > 250.0 ){
			auto f = fpaddles[ichan];
			f->SetParameter(0, max);
			f->SetParameter(1, xmax);
			f->SetParameter(2, 2.0);
			f->SetParameter(3, 100.0);
			h->Fit(f, "0Q");
			
			thit[ichan] = 4.0*f->GetParameter(1);
		}
	}
	
	for(int i=0; i<4; i++) if(thit[i]>0.0) htpaddles[i]->Fill(thit[i]);
	
	double tup = -1000.0;
	double tdown = -1000.0;
	double tmean_up = -1000;
	double tmean_dn = -1000;
	if( (thit[1]>0.0) && (thit[3]>0.0) ){
		tup = thit[1] - thit[3];
		tdiff_upstream->Fill( tup );
		tmean_up = (thit[1] + thit[3])/2.0;
	}
	if( (thit[0]>0.0) && (thit[2]>0.0) ){
		tdown = thit[0] - thit[2];
		tdiff_downstream->Fill( tdown );
		tmean_dn = (thit[0] + thit[2])/2.0;
	}
	if( (tup>-1000) && (tdown>-1000) ){
		tdiff_up_vs_down->Fill(tdown, tup);
		tup_minus_tdown->Fill(tup - tdown);
		tmean_diffs->Fill( (tup+tdown)/2.0 );
		
		for(auto iwire : fmwpchits) paddletime_vs_wire->Fill(iwire, (tup+tdown)/2.0);
	}
	if( (tmean_up>-1000) && (tmean_dn>-1000) ){
		tdiff_means->Fill( tmean_up - tmean_dn);
	}
	
	japp->RootFillUnLock(this);

	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t JEventProcessor_fmwpcana::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_fmwpcana::fini(void)
{
	for(auto h : hfmwpc ) delete h;
	for(auto h : hpaddles ) delete h;

	return NOERROR;
}