/**
* Rate
* author  Nagni
* version 1.0 - 27 April 2006
*
* Description: .
*
* Parameters:
* 	base - the path where to find the PAMELA unpacked root file.
*	outDirectory - the path where to put the output file.
*	orbPath - the path where to find an orb format for the output.
*
* version 1.0 - 27 April 2006
* First implementation
*/

#include <mcmd/McmdEvent.h>
#include <mcmd/McmdRecord.h>
#include <EventHeader.h>
#include <PscuHeader.h>
#include <TTree.h>
#include "cTle.h"
#include "cEci.h"
#include "cOrbit.h"
#include "TH2F.h"
#include "TFrame.h"
#include "TGraph.h"
#include "TCanvas.h"
#include "TASImage.h"
#include "TMarker.h"

#include "TString.h"
#include "TObjString.h"
#include "TStyle.h"
#include "TPaletteAxis.h"
#include "TROOT.h"
#include <sys/stat.h>

#define TRUE  1
#define FALSE 0
/**
 * 
 * @param base 
 * @param outDirectory 
 * @param xslPath 
 */
using namespace std;

void InitStyle() {
  gROOT->SetStyle("Plain");

  TStyle *myStyle[2], *tempo;
  myStyle[0]=new TStyle("StyleWhite", "white");
  myStyle[1]=new TStyle("StyleBlack", "black");

  tempo=gStyle;
  Int_t linecol, bkgndcol, histcol;
 
  for(Int_t style=0; style<2; style++) {

    linecol=kWhite*style+kBlack*(1-style);
    bkgndcol=kBlack*style+kWhite*(1-style);
    histcol=kYellow*style+kBlack*(1-style); // was 95

    myStyle[style]->Copy(*tempo);

    myStyle[style]->SetCanvasBorderMode(0);
    myStyle[style]->SetCanvasBorderSize(1);
    myStyle[style]->SetFrameBorderSize(1);
    myStyle[style]->SetFrameBorderMode(0);
    myStyle[style]->SetPadBorderSize(1);
    myStyle[style]->SetStatBorderSize(1);
    myStyle[style]->SetTitleBorderSize(1);
    myStyle[style]->SetPadBorderMode(0);
    myStyle[style]->SetPalette(1,0);
    myStyle[style]->SetPaperSize(20,27);	
    myStyle[style]->SetFuncColor(kRed);
    myStyle[style]->SetFuncWidth(1);
    myStyle[style]->SetLineScalePS(1); 
    myStyle[style]->SetCanvasColor(bkgndcol);
    myStyle[style]->SetAxisColor(linecol,"XYZ");
    myStyle[style]->SetFrameFillColor(bkgndcol);
    myStyle[style]->SetFrameLineColor(linecol);
    myStyle[style]->SetLabelColor(linecol,"XYZ");
    myStyle[style]->SetPadColor(bkgndcol);
    myStyle[style]->SetStatColor(bkgndcol);
    myStyle[style]->SetStatTextColor(linecol);
    myStyle[style]->SetTitleColor(linecol,"XYZ");
    myStyle[style]->SetTitleFillColor(bkgndcol);
    myStyle[style]->SetTitleTextColor(linecol);
    myStyle[style]->SetLineColor(linecol);
    myStyle[style]->SetMarkerColor(histcol);
    myStyle[style]->SetTextColor(linecol);
    
    myStyle[style]->SetGridColor((style)?13:kBlack);
    myStyle[style]->SetHistFillStyle(1001*(1-style));
    myStyle[style]->SetHistLineColor(histcol);
    myStyle[style]->SetHistFillColor((style)?bkgndcol:kYellow);
  }

  myStyle[1]->cd();

  gROOT->ForceStyle();

}

void Rate(TString base, TString outDirectory = "", TString format = "jpg", TString mapFile = ""){
	TTree                  *tr         = 0;
	pamela::McmdEvent      *mcmdev     = 0;
	pamela::McmdRecord     *mcmdrc     = 0;
	TArrayC                *mcmddata   = 0;
	ULong64_t        	nevents    = 0;
    	ULong64_t        	timesync   = 0;
    	pamela::EventHeader    *eh	   = 0;
  	pamela::PscuHeader     *ph 	   = 0;
  	stringstream 	oss;
    	double offsetTime = 0;
    	double absTime;
    	UInt_t i = 0;
    	UInt_t j = 0;
    	struct stat buf;
    	
    	i = stat (mapFile.Data(), &buf );
	// If the file does not exists
	if (i != 0){
		printf("The %s file does not exists.", mapFile.Data());
		exit(0);
	}
   	
    	TString filename = ((TObjString*)base.Tokenize('/')->Last())->GetString();
  	filename = ((TObjString*)filename.Tokenize('.')->First())->GetString();
    	// Test SGP4
   	string str1 = "SGP4 Test";
   	string str2 = "1 25544U 98067A   06117.32388940  .00009459  00000-0  64427-4 0  8131";
   	string str3 = "2 25544  51.6388  89.2928 0009043 155.3293 354.6512 15.75673697425143";

   	cTle tle1(str1, str2, str3);
   	cOrbit       orbit(tle1);
   	cEci         eci;
   	cCoordGeo    coo;
   	TH2F *rate   	= new TH2F("rate", base, 360, -180, 180, 161, -80.5, 80.5);
	TFile *rootFile = new TFile(base);
	if (rootFile->IsZombie()) {
		printf("The %s file does not exist", base.Data());
		exit(0);
	}
	
	
	/* 
	 * process Mcmd
	 */
	long int recEntries;
	tr = (TTree*)rootFile->Get("Mcmd");
	tr->SetBranchAddress("Mcmd", &mcmdev);
	nevents = tr->GetEntries();
	
	bool timeFound = FALSE;
	while ((i < nevents) | !timeFound){
		tr->GetEntry(i);
		recEntries = mcmdev->Records->GetEntries();
		while ((j < recEntries) | !timeFound){
			mcmdrc = (pamela::McmdRecord*)mcmdev->Records->At(j);
			if (mcmdrc->ID1 == 0xE0){
			mcmddata = mcmdrc->McmdData;
			timesync = (((ULong64_t)mcmddata->At(0)<<24)&0xFF000000) + 
					(((ULong64_t)mcmddata->At(1)<<16)&0x00FF0000) + 
					(((ULong64_t)mcmddata->At(2)<<8)&0x0000FF00)  + 
					(((ULong64_t)mcmddata->At(3))&0x000000FF);
			offsetTime = timesync - (mcmdrc->MCMD_RECORD_OBT)*(1./1000.);
			timeFound = TRUE;
			}
		j++;
		}
	i++;
	}
	
	tr = (TTree*)rootFile->Get("Physics");
	TBranch *headBr = tr->GetBranch("Header");
	tr->SetBranchAddress("Header", &eh);
	nevents = tr->GetEntries();
	//Fill variables from root-ple
  	for (i = 0; i < nevents; i++){
    		tr->GetEntry(i);
    		ph = eh->GetPscuHeader();
    		absTime = ((ph->GetOrbitalTime()*(1./1000.)) + absTime)/60;
    		orbit.getPosition(absTime, &eci);
      		coo = eci.toGeo();
      		rate->Fill(rad2deg(coo.m_Lon), rad2deg(coo.m_Lat));
      		/*
      		printf("             %16.8f %16.8f %16.8f\n",
             	rad2deg(coo.m_Lat),
             	rad2deg(coo.m_Lon),
             	coo.m_Alt);	
             	*/
    	}
	double posx=-1000,posy=-1000,oldposx=-1000,oldposy=-1000;

	int ptcnt=0, color=0;
	TMarker *tma=NULL;
	TLine *tli=NULL;
	double step=0;


	TImage *tImage=TImage::Open(mapFile);
	int width=(int)(tImage->GetWidth()*0.80);
	int height=(int)(tImage->GetHeight()*0.80);
	InitStyle();
	TCanvas *c1 = new TCanvas("c1","rate/orbit",-width, height); // - : removes the menu bars
	
	TH1F *hframe=NULL;
	hframe=gPad->DrawFrame(-180,-90,180,90);
	oss.str("");
	oss << filename << " - Event Rate (Hz)";
	hframe->SetTitle(oss.str().c_str());
	hframe->SetXTitle("Longitude (deg)");
	hframe->SetYTitle("Latitude (deg)");
	c1->cd();
	
	
	
	TPad *p2 = new TPad("p2", "p2", 0.10, 0.04, 0.983, 1);
	p2->Draw();
	p2->cd();
	TPaletteAxis *tpa=NULL;
	TH2F *forpal=new TH2F("forpal","",2,0,2,2,0,2);
	forpal->SetAxisColor(kBlack); //Delete the stat box
	forpal->SetStats(0); //Delete the stat box
	forpal->SetMinimum(0.1);
	forpal->SetMaximum(200);
	forpal->SetBinContent(5,1);    // just to initialize the histo
	forpal->SetContour(50);
	TPaletteAxis *tpp=(TPaletteAxis*)((forpal->GetListOfFunctions())->FindObject("palette"));
	step=forpal->GetContourLevel(1)-forpal->GetContourLevel(0);
	tpa=new TPaletteAxis(184,-90,195,90,forpal);
	tpa->SetLabelColor(kWhite);
	forpal->Draw("colz");

	c1->cd();
	TPad *p1 = new TPad("p1", "p1", 0.10,0.10,0.90,0.92);
	p1->Draw();
	p1->cd();
	tImage->Draw("X");
	c1->cd();
	gPad->RedrawAxis();
	c1->Update();
	
  c1->cd();
  ptcnt=0;
  tma=new TMarker();
  tma->SetMarkerStyle(4);
  tli=new TLine();
  tli->SetLineColor(kMagenta);
  Stat_t freq;
      for (Int_t kk = 0; kk < 360; kk++){
        for (Int_t jj = 0; jj < 161; jj++){
		freq = rate->GetBinContent(kk, jj);
		if (freq>0) {
			posx=(kk - 180); posy= jj - 80.5;
			// color: palette colors from 51 to 100 ie 50 levels
			color=51+(int) ((log10((rate==0)?0.1:freq)-log10(0.1))/step); // step defined by palette
			if (color>100) color=100; // just in case if max rate is not max...
			tma->SetMarkerColor(color);
			if (!(posx<0 && oldposx>0) && oldposx!=-1000 && oldposy!=-1000) tli->DrawLine(oldposx,oldposy,posx,posy); 
			tma->DrawMarker(posx,posy);
			oldposx=posx;
			oldposy=posy;	
		}
      	}
      }
  	oss.str("");
  	if (outDirectory == "") {
  		oss << filename.Data() << "_OrbitRate." << format.Data();
  	} else {
  		oss << outDirectory.Data() << filename.Data() << "_OrbitRate." << format.Data();
  	}
	c1->SaveAs(oss.str().c_str());
        rootFile->Close();
}

int main(int argc, char* argv[]){
    TString outDir  = "./";
    TString mapFile = "";
    TString format  = "jpg";
    
 if (argc < 2){
    printf("You have to insert at least the file to analyze and the mapFile \n");
    printf("Try '--help' for more information. \n");
    exit(1);
  }  

  if (!strcmp(argv[1], "--help")){
        printf( "Usage: OrbitRate FILE -map mapFile [OPTION] \n");
        printf( "mapFile have to be a mercator map image [gif|jpg|png] \n");
        printf( "\t --help                 Print this help and exit \n");
        printf( "\t -outDir[path]          Path where to put the output [default ~/tmp] \n");
        printf( "\t -format[jpg|gif|ps]    Format for output files [default 'jpg'] \n");
        exit(1);
  }

  for (int i = 2; i < argc; i++){
    if (!strcmp(argv[i], "-outDir")){
        if (++i >= argc){
            printf( "-outDir needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
            outDir = argv[i];
            continue;
        }
    }

    if (!strcmp(argv[i], "-map")){
        if (++i >= argc){
            printf( "-map needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
            mapFile = argv[i];
            continue;
        }
    }

    if (!strcmp(argv[i], "-format")) {
    	if (++i >= argc){
            printf( "-format needs arguments. \n");
            printf( "Try '--help' for more information. \n");
            exit(1);
        } else {
        	format = argv[i];
        	continue;
        }
    }
}
if (mapFile != ""){
	Rate(argv[1], outDir, format, mapFile);
} else {
	printf("You have to insert at least the file to analyze and the mapFile \n");
	printf("Try '--help' for more information. \n");
    }
}