//********************************
//
// Class CalFactFits
// Wrapper class for factfits.h (or fits.h, pyfits.h)
// Provides fast access to calibrated events of FACT raw files
//
// written by Thomas Kraehenbuehl, ETH Zurich
// tpk@phys.ethz.ch, +41 44 633 3973
// April 2012
//
//********************************
//
// Compilation (root in pyfact directory)
// root [0] gSystem->Load("/usr/lib64/libz.so");
// root [1] .L izstream.h++
// root [2] .L factfits.h++
// root [3] .L calfactfits.h++
// 
// Try /usr/lib/x86_64-linux-gnu/libz.so at ETH for the libz.so file.
//
// Note: the produced izstream_h.so must be in LD_LIBRARY_PATH, e.g.
// export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/home/guest/kraehenb/software/PythonIncludes/
//
// Usage in Python:
// See pyscripts/examples/CalFitsTest.py
//
//********************************

//ToDo: shared library creation debuggen

#ifndef CALFACTFITS_H
#define CALFACTFITS_H

#include <cstdio>
#include <string>

#ifndef __MARS__
#include <vector>
#include <iomanip>
#include <iostream>
#define gLog cerr
#define ___err___ ""
#define ___all___ ""
#else
#include "MLog.h"
#include "MLogManip.h"
#define ___err___ err
#define ___all___ all
#endif

#ifdef __EXCEPTIONS
#include <stdexcept>
#endif

#include "extern_Mars_mcore/factfits.h"

class CalFactFits
{
public:
	//No standard constructor CalFactFits()!
	
	//Direct handlers of the files
	factfits datafile, calibfile; //Class name should be PyFits or better FactPyFits...
	
	//Basic file parameters
	UInt_t calib_nroi, calib_npix;
	UInt_t calib_blm_size, calib_gm_size, calib_tom_size;
	UInt_t data_nroi, data_npix, data_ndata;
	
	//Common variables
	UInt_t nroi, npix, nevents;
	
	//Calibration variables
	float* calib_baselinemean;
	float* calib_gainmean;
	float* calib_triggeroffsetmean;
    // The same variables as above, but this time
    // they are doubled in memory, so the modulo operation is not necessary.
    // idea by TPK again
    // since we like to work with doubles for some reason, I don't know yet
    // I cast the float to double already in this stage.
    
    double* baseline;
    double* gain;
    double* trigger_offset;
	//Using <vector> instead of arrays makes no visible difference
	//ToDo: use arrays of size 1440x1024 (x2 for wrap-arounds) and read all variables into those
	
	//Event variables
	UInt_t event_id;
	UShort_t event_triggertype;
	int16_t* event_data;
	int16_t* event_offset;
	int32_t* event_boardtimes;
	double* npcaldata;
	
	CalFactFits(const string &datafilename, const string &calibfilename) //Constructor with two filenames
		: datafile(datafilename),
			calibfile(calibfilename),
			npcaldata(NULL)
	{
        //cout << "Constructor called" << endl;
		//Read basic parameters of the two files
//		std::cout << "...Reading basic file parameters..." << std::endl;
		try {
			calib_nroi = calibfile.GetUInt("NROI");
		}
		catch(runtime_error &err) {
			std::cout << err.what() << std::endl;
			std::cout << "Standard value NROI=1024 assumed." << std::endl;
			calib_nroi = 1024;
		}
		try {
			calib_npix = calibfile.GetUInt("NPIX");
		}
		catch(runtime_error &err) {
			std::cout << err.what() << std::endl;
			std::cout << "Standard value NPIX=1440 assumed." << std::endl;
			calib_npix = 1440;
		}
		
		data_nroi = datafile.GetUInt("NROI");
		data_npix = datafile.GetUInt("NPIX");
		data_ndata = datafile.GetN("Data");
		
		calib_blm_size = calibfile.GetN("BaselineMean")/calib_npix;
		calib_gm_size = calibfile.GetN("GainMean")/calib_npix;
		calib_tom_size = calibfile.GetN("TriggerOffsetMean")/calib_npix;
		
//		std::cout << "Column sizes: " << calib_blm_size << " " << calib_gm_size << " " << calib_tom_size << std::endl;
		
		//Define the common variables
		if((calib_nroi==data_nroi)&&(calib_npix==data_npix)&&(data_nroi*data_npix==data_ndata)&&(calib_blm_size==calib_gm_size)&&(calib_tom_size==calib_nroi)) {
			nroi = data_nroi;
			npix = data_npix;
		}
		else {
			ostringstream str;
			str << "Data/calib file error: NROI mismatch, NPIX mismatch, data column size wrong or calib columns mismatch.";
#ifdef __EXCEPTIONS
			throw runtime_error(str.str());
#else
			gLog << ___err___ << "ERROR - " << str.str() << endl;
			return;
#endif
		}
		nevents = datafile.GetNumRows();
		
		//Read the calibration data
//		std::cout << "...Reading calibration data..." << std::endl;
		calib_baselinemean = new float[calibfile.GetN("BaselineMean")];
        calibfile.SetPtrAddress("BaselineMean", calib_baselinemean, calibfile.GetN("BaselineMean"));
        baseline = new double[calibfile.GetN("BaselineMean")*2];
        
		calib_gainmean = new float[calibfile.GetN("GainMean")];
		calibfile.SetPtrAddress("GainMean", calib_gainmean, calibfile.GetN("GainMean"));
        gain = new double[calibfile.GetN("GainMean")*2];
        
		calib_triggeroffsetmean = new float[calibfile.GetN("TriggerOffsetMean")];
		calibfile.SetPtrAddress("TriggerOffsetMean", calib_triggeroffsetmean, calibfile.GetN("TriggerOffsetMean"));
		trigger_offset = new double[calibfile.GetN("TriggerOffsetMean")*2];
        
        calibfile.GetRow(0);
        
int orig_index, new_index1, new_index2;
        for (int pix=0 ; pix < (int)calib_npix; ++pix)
        {
            for (int sl = 0; sl < (int)calib_blm_size; ++sl)
            {
                orig_index = pix * (int)calib_blm_size + sl;
                new_index1 = 2*pix * (int)calib_blm_size + sl;
                new_index2 = (2*pix+1) * (int)calib_blm_size + sl;
                
                baseline[new_index2] = baseline[new_index1] = double(calib_baselinemean[orig_index]) *4096./2000.;
            }
        }
        for (int pix=0 ; pix < (int)calib_npix; ++pix)
        {
            for (int sl = 0; sl < (int)calib_gm_size; ++sl)
            {
                orig_index = pix * (int)calib_gm_size + sl;
                new_index1 = 2*pix * (int)calib_gm_size + sl;
                new_index2 = (2*pix+1) * (int)calib_gm_size + sl;
                
                gain[new_index2] = gain[new_index1] = double(calib_gainmean[orig_index]) /( 1907.35 /4096. *2000.);
            }
        }
        
        for (int pix=0 ; pix < (int)calib_npix; ++pix)
        {
            for (int sl = 0; sl < (int)calib_tom_size; ++sl)
            {
                orig_index = pix * (int)calib_tom_size + sl;
                new_index1 = 2*pix * (int)calib_tom_size + sl;
                new_index2 = (2*pix+1) * (int)calib_tom_size + sl;
                
                trigger_offset[new_index2] = trigger_offset[new_index1] = double(calib_triggeroffsetmean[orig_index]) *4096. /2000.;
            }
        }
        
        
		
		//Set the event pointers
//		std::cout << "...Setting event pointers..." << std::endl;
		datafile.SetRefAddress("EventNum", event_id);
		datafile.SetRefAddress("TriggerType", event_triggertype);
		
		event_data = new int16_t[data_ndata];
		datafile.SetPtrAddress("Data", event_data, data_ndata);
		
		event_offset = new int16_t[datafile.GetN("StartCellData")];
		datafile.SetPtrAddress("StartCellData", event_offset, datafile.GetN("StartCellData"));
		
		event_boardtimes = new int32_t[datafile.GetN("BoardTime")];
		datafile.SetPtrAddress("BoardTime", event_boardtimes, datafile.GetN("BoardTime"));
	}
	
	~CalFactFits() //Standard destructor
	{
		delete[] calib_baselinemean;
		delete[] calib_gainmean;
		delete[] calib_triggeroffsetmean;
		delete[] event_data;
		delete[] event_offset;
		delete[] event_boardtimes;
        delete[] baseline;
        delete[] gain;
        delete[] trigger_offset;
	}
	
	bool GetCalEvent() //Read calibrated event into the event variables
	{
		if(!npcaldata) {
			ostringstream str;
			str << "Pointer to the calibrated data not initialized!";
#ifdef __EXCEPTIONS
			throw runtime_error(str.str());
#else
			gLog << ___err___ << "ERROR - " << str.str() << endl;
			return false;
#endif
		}
		if(datafile.GetNextRow() == false) {
//			std::cout << "Last event reached..." << std::endl;
			return false;
		}
		else 
        {
			UInt_t drs_calib_offset;
            double raw, raw_bsl, shifted;
			for(UInt_t pixel = 0; pixel < data_npix; pixel++)
            {
				for(UInt_t slice = 0; slice < data_nroi; slice++)
                {
                    drs_calib_offset = slice + event_offset[pixel];
                    raw = (double)event_data[pixel*data_nroi+slice];
                    raw_bsl = raw - baseline[pixel*calib_blm_size*2+drs_calib_offset];
                    shifted = raw_bsl - trigger_offset[pixel*data_nroi*2+slice];
                    npcaldata[pixel*data_nroi+slice] = shifted / gain[pixel*calib_blm_size*2+drs_calib_offset];
					//Note: data_nroi=calib_nroi, calib_blm_size=calib_gm_size
				}

                // stolen from TBs MARS ... DrsCalib.h or something like that
                 double *p = npcaldata+(pixel*data_nroi);

                // second loop for despiking
                for (size_t i=1; i<data_nroi-2; i++)
                {
                    if (p[i]-p[i-1]>25 && p[i]-p[i+1]>25)
                    {
                        //std::cout << "single cand:\t" << pixel << "\t" << i << std::endl;
                        p[i] = (p[i-1]+p[i+1])/2;
                    }
                    if (p[i]-p[i-1]>22 && fabs(p[i]-p[i+1])<4 && p[i+1]-p[i+2]>22)
                    {
                        //std::cout << "double cand:\t" << pixel << "\t" << i << std::endl;
                        p[i] = (p[i-1]+p[i+2])/2;
                        p[i+1] = p[i];
                    }
                }
                
			}
		}
		return true;
	}
	
	void SetNpcaldataPtr(double *numpyptr) //Set the pointer for the calibrated data to the numpy array
	{
		npcaldata = numpyptr;
		return;
	}
};
#endif /* CALFACTFITS_H */