// $Id$
//
//    File: MyProcessor.cc
// Created: Mon Aug 21 14:42:12 EDT 2006
// Creator: davidl (on Darwin swire-b241.jlab.org 8.7.0 powerpc)
//

#include <fstream>
using namespace std;

#include <TF1.h>
#include <JANA/JApplication.h>

#include "DEventSourceET/DADC.h"
#include "DEventSourceET/DTrigger.h"

#include "MyProcessor.h"

//------------------
// init
//------------------
jerror_t MyProcessor::init(void)
{
	// pedestal histograms
	channel_vs_ped_slot18 = new TH2F("channel_vs_ped_slot18", "ADC pedestal vs. channel slot 18", 512, 0.0, 512.0, 32, -0.5, 31.5);
	channel_vs_ped_slot19 = new TH2F("channel_vs_ped_slot19", "ADC pedestal vs. channel slot 19", 512, 0.0, 512.0, 32, -0.5, 31.5);

	Nevents_processed=0;
	runnumber=0;

	return NOERROR;
}

//------------------
// brun
//------------------
jerror_t MyProcessor::brun(JEventLoop *eventLoop, int runnumber)
{
	if(runnumber>this->runnumber)this->runnumber=runnumber;

	return NOERROR;
}

//------------------
// evnt
//------------------
jerror_t MyProcessor::evnt(JEventLoop *loop, int eventnumber)
{
	// Get data
	vector<const DADC*> dadcs;
	vector<const DTrigger*> dtriggers;
	loop->Get(dadcs);
	loop->Get(dtriggers);

	// We need non-beam related triggers
	if(dtriggers.size()){
		int trig_mask = (0x1<<6)+(0x1<<5); // pulser or cosmic trigger
		if((dtriggers[0]->latch&trig_mask)==0)return NOERROR;
	}

	// Fill histograms
	for(unsigned int i=0; i<dadcs.size(); i++){
		const DADC *dadc = dadcs[i];
		
		if(dadc->overflow)continue;
		
		if(dadc->slot==18)channel_vs_ped_slot18->Fill(dadc->adc, dadc->channel);
		if(dadc->slot==19)channel_vs_ped_slot19->Fill(dadc->adc, dadc->channel);
		
	}
	
	// Only need first 1000 good events
	if(Nevents_processed++>=1000)loop->GetJApplication()->Quit();
	
	return NOERROR;
}

//------------------
// erun
//------------------
jerror_t MyProcessor::erun(void)
{	
	return NOERROR;
}

//------------------
// fini
//------------------
jerror_t MyProcessor::fini(void)
{
	cout<<" Found "<<Nevents_processed<<" pulser/cosmic events"<<endl;

	char fname[256];
	sprintf(fname, "ped_slot18_run%05d.txt",runnumber);
	WritePedFile(fname,channel_vs_ped_slot18);

	sprintf(fname, "ped_slot19_run%05d.txt",runnumber);
	WritePedFile(fname,channel_vs_ped_slot19);

	return NOERROR;
}

//------------------
// WritePedFile
//------------------
void MyProcessor::WritePedFile(const char* fname, TH2F *chan_vs_ped)
{
	// open output file
	ofstream fout(fname);
	
	// Loop over ADC channels, projecting them out and fitting them
	TAxis *xaxis = chan_vs_ped->GetXaxis();
	TH1D *ped = new TH1D("ped","Pedestal", xaxis->GetNbins(), xaxis->GetXmin(), xaxis->GetXmax());
	for(int i=1; i<=32; i++){
		chan_vs_ped->ProjectionX("ped", i,i);
		float max = ped->GetBinCenter(ped->GetMaximumBin());
		ped->Fit("gaus","Q","",max-5, max+5);
		TF1 *func = ped->GetFunction("gaus");
		float amp = func->GetParameter(0);
		float mean = func->GetParameter(1);
		float sigma = func->GetParameter(2);
		fout<<(int)mean<<endl;
		cout<<"channel "<<i-1<<":  amp:"<<amp<<"  mean:"<<mean<<"   sigma:"<<sigma<<endl;
	}
	delete ped;

	// close output file
	fout.close();
}