/*
GetDrsOffset.C - a ROOT macro calculating the DRS Offset calibration constants

    -- TESTING --

 In case you have taken a file, which is good for calibrating the DRS Offset,
  then you can calculate the offset using this macro.
  It just calculates the mean value and the RMS for all DRS bins of all pixels
  in the file, for the number of events you like.

  The output is two TH2Fs, which contain all requested data on the one hand and
  on the other hand are good for having an overview.

  call it like this:
  GetDrsOffset("path/20111111_001.fits.gz")

  GetDrsOffset("path/20111111_001.fits.gz", "foo/output.root")
  saves the output to a root file

  GetDrsOffset("path/20111111_001.fits.gz", "", "bar/test.bin")
  saves the output to a binary file.

  GetDrsOffset("path/20111111_001.fits.gz", "foo/out.root", "bar/test.bin", false, -1)
  produces all output, but does not open a ROOT Canvas.
  -1 : means process all events in the given file.

*/


#include <TROOT.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <Getline.h>
#include <TFile.h>
#include <TStyle.h>

#include <stdio.h>
#include <stdint.h>
#include <cstdio>
#include <fstream>
#include <string.h>

#define HAVE_ZLIB
#include "fits.h"

#include "openFits.h"
#include "openFits.c"


int NEvents;
vector<int16_t> AllPixelDataVector;
vector<int16_t> StartCellVector;
unsigned int CurrentEventID;
size_t PXLxROI;
UInt_t RegionOfInterest;
UInt_t NumberOfPixels;

TH2F * hOffset;
TH2F * hOffsetRMS;
double mean[1440][1024]={{0}};
int n[1440][1024]={{0}};
double rms[1440][1024]={{0}};
TObjArray hList;


int GetDrsOffset(
	char *datafilename 		= "data/20111016_013.fits.gz",
	const char *RootFileName = "",
	const char *OutFileName = "",
	bool graphic = true,
	int nevents 			= -1,
	int	verbosityLevel = 1
	)
{
gStyle->SetPalette(1,0);


	fits * datafile;
	NEvents = openDataFits( datafilename, &datafile,
		AllPixelDataVector, StartCellVector, CurrentEventID,
		RegionOfInterest, NumberOfPixels, PXLxROI, verbosityLevel);
	if (NEvents == 0){
		cout << "return code of OpenDataFile:" << datafilename<< endl;
		cout << "is zero -> aborting." << endl;
		return 1;
	}

	if (verbosityLevel > 0)
		cout <<"number of events in file: "<< NEvents << "\t";
	if ( nevents == -1 || nevents > NEvents ) nevents = NEvents; // -1 means all!
	if (verbosityLevel > 0)
		cout <<"of, which "<< nevents << " will be processed"<< endl;

	int npixel = -1;
	if (verbosityLevel > 0)
		cout <<"Total # of Pixel: "<< NumberOfPixels << "\t";
	if ( npixel == -1 || npixel > (int)NumberOfPixels  ) npixel = (int)NumberOfPixels; // -1 means all!
	if (verbosityLevel > 0)
		cout <<"of, which "<< npixel << " will be processed"<< endl;

	
	if (verbosityLevel > 0)
		cout <<" The ROI of this datafile is:" << RegionOfInterest << endl;

	double cur;
  // -------------------- OFFSET --------------------------------------------
  for ( int ev = 0; ev < nevents; ev++) {
	printf("%5d / %5d\b\b\b\b\b\b\b\b\b\b\b\b\b",ev,nevents);
    datafile->GetRow( ev );
    for ( int pix = 0 ; pix < npixel; pix++ ){
      for (unsigned int sl=0; sl< RegionOfInterest; sl++){
        cur = AllPixelDataVector[pix * RegionOfInterest + sl];
        mean[pix][(sl+StartCellVector[pix])%1024] += cur;
				++n[pix][(sl+StartCellVector[pix])%1024];
        rms[pix][(sl+StartCellVector[pix])%1024] += cur * cur;
      }
    }
  }
  for (int p=0; p<1440; p++)
    for (int s=0; s<1024; s++){
			if (n[p][s]==0){
				cout << "Error: no data for pixel:" << p <<" slice:" << s << endl;
				return 1;
			}
  	  mean[p][s]/=n[p][s];
   	 rms[p][s]/=n[p][s];
   	 rms[p][s]-=mean[p][s]*mean[p][s];
   	 rms[p][s]=sqrt(rms[p][s]);
  	}


	// ----------------- For plotting only -------------------
	if (graphic){

	TCanvas *cOffset = NULL;
	cOffset = new TCanvas("cOffset","Offset",10,10,800,400);
	cOffset->Divide(1, 2);

  hOffset = new TH2F("hOffset","Offset",1440,-0.5,1439.5, 1024, -0.5, 1023.5);
  hOffset->GetXaxis()->SetTitle( "pixel" );
  hOffset->GetYaxis()->SetTitle( "slice" );
  hList.Add( hOffset );

  hOffsetRMS = new TH2F("hOffsetRMS","Offset RMS",1440,-0.5,1439.5, 1024, -0.5, 1023.5);
  hOffsetRMS->GetXaxis()->SetTitle( "pixel" );
  hOffsetRMS->GetYaxis()->SetTitle( "slice" );
  hList.Add( hOffsetRMS );

    for ( int pix = 0 ; pix < 1440; pix++ ){
      for (unsigned int sl=0; sl< 1024; sl++){
        hOffset->Fill(pix,sl,mean[pix][sl]);
        hOffsetRMS->Fill(pix,sl,rms[pix][sl]);
      }
    }
    cOffset->cd(1);
    hOffset->Draw("COLZ");
    cOffset->cd(2);
    hOffsetRMS->Draw("COLZ");
    cOffset->Modified();
    cOffset->Update();
	}

	if ( strlen(OutFileName) > 0 ) {
		ofstream out(OutFileName, ios::out | ios::binary);
  	if(!out) {
  		cout << "Cannot open file." << OutFileName << endl;
   	 return 1;
  	}
   	  out.write((char *) &mean, sizeof mean);
   	  out.write((char *) &rms, sizeof rms);
   	  out.close();
		cout << "binary output written to \"" << OutFileName << "\"" << endl;
	}else
		cout << "no output generated" << endl;
   
/*
  double offset[1440][1024];
  ifstream in(OutFileName, ios::in | ios::binary);
  in.read((char*) &offset, sizeof offset);
  in.close();
	
	
    for ( int pix = 0 ; pix < 1440; pix++ ){
      for (unsigned int sl=0; sl< 1024; sl++){
				cout << offset[pix][sl] << endl;
      }
    }
*/
	
	TFile tf( RootFileName, "RECREATE");
	hList.Write(); // write the major histograms into the top level directory
	tf.Close(); // close the file
	return( 0 );
}