/********************************************************************
 *
 * Create anisotropy search skymap, 
 *
 * Adopted some Auger script for Halo anaysis 
 * Carsten Rott (CCAPP) 2010
 *  
 *******************************************************************/


#include <iostream>
#include <iomanip>

#include <stdlib.h>
#include <string>
#include <stdio.h>
#include <ctype.h>
#include <math.h>
#include "TGaxis.h"
  

#include "TLegend.h"
#include "TBasket.h" 
#include "Riostream.h"
#include "TEventList.h"
#include "TFile.h"
#include "TSystem.h"
#include "TChain.h"
#include "TTree.h"
#include "TROOT.h"
#include "TLatex.h"


#include "TMath.h"
#include "TObject.h"
#include "TObjString.h"
#include "TObjArray.h"
#include "TTimeStamp.h"
#include "TFile.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TF1.h"
#include "TStyle.h"
#include "TLegend.h"
#include "TNamed.h"

#include "plot_tools.h"

/*
  root -l
  gSystem->CompileMacro("SkymapCode_DM_on_off_source_forPRD.C");
  SkymapCode_DM_on-off-source_forPRD()

  or alternatively 
  root -l SkymapCode_DM_on_off_source_forPRD.C

 */


void SkymapCode_DM_on_off_source_forPRD() {

    plot_nice();

    //gStyle->SetPalette(2);
    //gStyle->SetNumberContours(77);
    set_my_palette(3); // from plot_tools.h


  //example histogram
//  TFile *file100 = TFile::Open("SkyMap_output.root");
  TFile *file100 = TFile::Open("SkyMap_HaloFlux_output_.root");
  TH2F* flux = (TH2F*)file100->Get("h2_skymap_theta_RA_signal");
//  TFile *file100 = TFile::Open("exampleFile.root");
//  TH2D* flux     = (TH2D*)file100->Get("flux_cr11");

  //flux->Scale(206);

  float pi=3.1415926; 
  float conv=pi/180; // I am also aware of TMath::DegToRad() and TMath::Pi() which could be used there...
  

  /// Draw the grid ///

  const int Nl=7; // Number of drawn latitudes
  const int NL=13; // Number of drawn longitudes
  const int NC=8; // Number of circles drawn around the galactic center 
  int M=60; // number of points in a graph
  int M2=3600; // number of points in a graph

  TGraph   *latitudes[Nl];
  TGraph   *longitudes[NL];
  TGraph   *circles[NC];
  TGraph   *circles_dotted[NC];

  TGraph   *circles_off[NC];
  TGraph   *circles_dotted_off[NC];
  
  for(int j=0;j<Nl;++j)
    {
      latitudes[j]=new TGraph();
      float la=-90+180/(Nl-1)*j;
      
      for(int i=0;i<M+1;++i)
	{
	  float lo= -180+360/M*i;
	  float z = sqrt(1+cos(la*conv)*cos(lo*conv/2));
	  float x = 180*cos(la*conv)*sin(lo*conv/2)/z;
	  float y = 90*sin(la*conv)/z;
	  latitudes[j]->SetPoint(i,x,y);
	}
    } 
  
  for(int j=0;j<NL;++j)
    {
      longitudes[j]=new TGraph();
      float lo=-180+360/(NL-1)*j;
      
      for(int i=0;i<M+1;++i)
	{
	  float la= -90+180/M*i;
	  float z = sqrt(1+cos(la*conv)*cos(lo*conv/2));
	  float x = 180*cos(la*conv)*sin(lo*conv/2)/z;
	  float y = 90*sin(la*conv)/z;
	  longitudes[j]->SetPoint(i,x,y);
	}
    }
 

  /// 10 deg circles around galactic center on Northern Hemisphere

  for(int j=0;j<NC;++j) // j=0,1,2, ... =0,10,20,... degrees (10 degree steps)
    {
      circles[j]=new TGraph();
      circles_dotted[j]=new TGraph();
      circles_off[j]=new TGraph();
      circles_dotted_off[j]=new TGraph();
      //float lo=-180+360/(NL-1)*j;
      //float lo=(j+1)*10-180+244;
      
      int horizon_count =0; // number of "grid" points above the horizon ... used for circles/circles_dotted

      for(int i=0;i<M2+1;++i)
	{
	    float lo= -30+180-244+((j+1)*10)*cos(float(i)/M2*pi*2);  // [-180,180]
	    //lo=-lo; // invert it
	    float la= -28+((j+1)*10)*sin(float(i)/M2*pi*2);        //   [-90,90]
	    float z = sqrt(1+cos(la*conv)*cos(lo*conv/2));
	    float x = 180*cos(la*conv)*sin(lo*conv/2)/z;
	    float y = 90*sin(la*conv)/z;

	    float z_off = sqrt(1+cos(la*conv)*cos((lo+180)*conv/2));
	    float x_off = 180*cos(la*conv)*sin((lo+180)*conv/2)/z_off;
	    float y_off = 90*sin(la*conv)/z_off;
	    if(la>-6.5) // 5 degree below the horizon is still drawn
	    {
		circles[j]->SetPoint(i-horizon_count,x,y); // if about the horizon
		circles_off[j]->SetPoint(i-horizon_count,x_off,y_off); // if about the horizon
		//printf("[j=%d] = (i-hor)=%d, i=%d\n",j,i-horizon_count,i);
	    }
	    else
	    {		
		circles_dotted[j]->SetPoint(horizon_count,x,y);
		circles_dotted_off[j]->SetPoint(horizon_count,x_off,y_off);
		horizon_count++;
	    }
	}
    }

  // labels and text //

  gStyle->SetPadBottomMargin(.08);
  gStyle->SetPadTopMargin(0.08);
  gStyle->SetPadLeftMargin(0.04);
  gStyle->SetPadRightMargin(0.15);
  Double_t ymin = -90;
  Double_t ymax = 90;
  Double_t dy = (ymax-ymin)/0.8; //10 per cent margins top and bottom
  Double_t xmin = -180;
  Double_t xmax = 180;
  Double_t dx = (xmax-xmin)/0.8; //10 per cent margins left and right
  //double conts2[] = {0.001,0.25,0.5,0.75,1.,1.25,1.5,1.75,2.,2.25,2.5,2.75,3.,3.25,3.5,3.75,4.,4.25,4.5,4.75,5.,5.25, 5.5, 5.75 };
  double x1[] = {0.02,-0.29,-0.51,-0.51,-0.28,0.02,-0.6};
  double y1[] = {-.61,-.53,-.3,.3,.52,.6,0.03};
  TLatex lat;
  lat.SetTextSize(0.05);
  lat.SetTextFont(42);
  lat.SetTextAlign(32);  

  TCanvas *smoothcan21 = new TCanvas("smoothcan21","Draw the skymap ");
  smoothcan21->SetFillColor(0);
  smoothcan21->SetBorderSize(0);
  smoothcan21->SetFrameBorderMode(0);
  smoothcan21->SetFrameLineColor(2); 
  smoothcan21->SetFrameBorderMode(0);    
  

  //set the contour levels
  //flux->SetContour(22,conts2);
  //draw with ROOT's aitoff projection (other projections are available, see ftp://root.cern.ch/root/doc/3Histograms.pdf, p31 (pdf doc page 9)
  flux->SetLineColor(3);
  flux->Draw("z aitoff");//zcol");
  flux->SetTitle("");
  
//  flux->GetZaxis()->SetRangeUser(0.,5.75);
  smoothcan21->Range(0.0, 0.0, 1.0, 1.0);
  flux->GetYaxis()->SetLabelOffset(999);
  flux->GetYaxis()->SetTickLength(0);
  flux->GetXaxis()->SetLabelOffset(999);
  flux->GetXaxis()->SetTickLength(0);
  gPad->Update();
 
  // labels //

  lat.DrawLatex(x1[0],y1[0]," -90#circ");
  lat.DrawLatex(x1[1],y1[1]," -60#circ");
  lat.DrawLatex(x1[2],y1[2]," -30#circ");
  lat.DrawLatex(x1[3],y1[3]," 30#circ");
  lat.DrawLatex(x1[4],y1[4]," 60#circ");
  lat.DrawLatex(x1[5],y1[5]," 90#circ");
  lat.DrawLatex(x1[6],y1[6]," 0#circ");

  TPad *pad21 = new TPad("pad21","",0,0,1,1);
  pad21->SetFillStyle(4000);
  pad21->SetFillColor(0);
  pad21->SetBorderSize(0);
  pad21->SetFrameBorderMode(0);
  pad21->SetFrameLineColor(0); 
  pad21->SetFrameBorderSize(0);
  pad21->Draw();
  pad21->cd();
 
  //pad offsets are set by .rootstyle.C file; all values are +0.01
  pad21->Range(xmin-0.04*dx,ymin-0.08*dy,xmax+0.15*dx,ymax+0.08*dy);
  for(int j=0;j<Nl;++j)
  {      
      latitudes[j]->SetLineColor(1);// black = 1 ; white = 0;
      latitudes[j]->SetLineStyle(3);
      latitudes[j]->Draw("c");
  }
  for(int j=0;j<NL;++j) 
  {
      longitudes[j]->SetLineColor(1); // black = 1 ; white = 0;
      longitudes[j]->SetLineStyle(3);
      longitudes[j]->Draw("c"); 
  }
  /// done with the grid ... now do the circle

    int circles_line_color=12; // black = 1 ; white = 0; red = 2

  for(int j=0;j<NC;++j)  
  {
      circles[j]->SetLineColor(circles_line_color); 
      circles[j]->SetLineStyle(2); 
      if(j==NC-1) circles[j]->SetLineStyle(1); // solid line
      circles[j]->Draw("c"); 

      circles_off[j]->SetLineColor(circles_line_color); 
      circles_off[j]->SetLineStyle(2); 
      if(j==NC-1) 
      {
	  circles_off[j]->SetLineStyle(1); // solid line
	  circles_off[j]->Draw("c"); 
      }
  }
  //for(int j=0;j<NC;++j)
  for(int j=0;j<4;++j)
  {
      circles_dotted[j]->SetLineColor(circles_line_color); 
      circles_dotted[j]->SetLineStyle(4);
      circles_dotted[j]->Draw("l"); 
  }     

  // now add more labels:
  TLatex lat2;
  lat2.SetTextSize(0.05);
  lat2.SetTextFont(42);
  lat2.SetTextAlign(32);  
  lat2.DrawLatex(-.3,-.05," #Delta RA=0#circ");  
  lat2.DrawLatex(.3,-.05," #Delta RA=180#circ");  
  lat2.DrawLatex(-.3,-.2," GC");
  TLatex lat3;
  lat3.SetTextSize(0.05);
  lat3.SetTextFont(42);
  lat3.SetTextAlign(32);    
  lat3.SetTextColor(0);    
  lat3.DrawLatex(-.3,.3," on source");  
  lat3.DrawLatex(0.3,.3," off source");  

}