#include <iostream>
#include <iomanip>

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

#include "TBasket.h" 
#include "Riostream.h"
#include "TSystem.h"
#include "TROOT.h"

#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TStyle.h"
#include "TNamed.h"




void SkymapProb() {
  const double Pi = 3.1415; 
  
  //effective area, efficiency histograms.  1 deg^2 binning
  TFile *file2 =  TFile::Open("SkymapProb.root "); 
  
  // area, efficiency for cr source flux
  TH2D *eff_dist_cel = (TH2D*)file2->Get("efficiency_dist_cel"); 
  //cr flux limit histogram
  TH2D* flux_limit_cr1  = new TH2D("flux_limit_cr1","cr flux limit",360,-180,180,180,-90.,90.); 

  for (int i = 1; i <= eff_dist_cel->GetNbinsX(); i++){
    for (int j = 1; j <= eff_dist_cel->GetNbinsY(); j++){
      double eff_val = eff_dist_cel->GetBinContent(i,j);
      if(eff_val>0.) flux_limit_cr1->SetBinContent(i,j,0.05/eff_val);    
      else flux_limit_cr1->SetBinContent(i,j,0.);    

    }
  }
  
  float conv=3.1415926/180; // I am also aware of TMath::DegToRad() and TMath::Pi() which could be used there...
  
  const int Nl=19; // Number of drawn latitudes
  const int NL=19; // Number of drawn longitudes
  int M=30;

  TGraph   *latitudes[Nl];
  TGraph   *longitudes[NL];

  TGraph *try1 = new TGraph();  
  TGraph *try2 = new TGraph();  

  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);
	  if(j==1)  try1->SetPoint(i,x,y);
	  if(j==5)  try2->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);
	}
    }

  TCanvas *CanCanS = new TCanvas("CanCanS","skymap without aitoff");
  flux_limit_cr1->Draw("zcol ");

  for(int j=0;j<Nl;++j) latitudes[j]->Draw("l same");
  for(int j=0;j<NL;++j) longitudes[j]->Draw("l same"); 

  TCanvas *smoothcan2 = new TCanvas("smoothcan2","skymap with aitoff");
  flux_limit_cr1->Draw("z aitoff  ");
  for(int j=0;j<Nl;++j) latitudes[j]->Draw("l same");
  for(int j=0;j<NL;++j) longitudes[j]->Draw("l same"); 
  
}