#include <TMath.h>

aitoff2xy(Double_t l, Double_t b, Double_t &Al, Double_t &Ab)
{
   Double_t x, y;

   Double_t alpha2 = (l/2)*TMath::DegToRad();
   Double_t delta  = b*TMath::DegToRad();
   Double_t r2     = TMath::Sqrt(2.);
   Double_t f      = 2*r2/TMath::Pi();
   Double_t cdec   = TMath::Cos(delta);
   Double_t denom  = TMath::Sqrt(1. + cdec*TMath::Cos(alpha2));
   if (denom == 0) return;
   if (f == 0) return;
   x      = cdec*TMath::Sin(alpha2)*2.*r2/denom;
   y      = TMath::Sin(delta)*r2/denom;
   x     *= TMath::RadToDeg()/f;
   y     *= TMath::RadToDeg()/f;
   //  x *= -1.; // for a skymap swap left<->right
   Al = x;
   Ab = y;
}

void SkymapProb() {
  
   //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=TMath::Pi()/180; // I am also aware of TMath::DegToRad() and TMath::Pi() which could be used there...
   float la, lo, x, y, z;
  
   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];

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

   TCanvas *c = new TCanvas("c","skymap with aitoff");
   flux_limit_cr1->Draw("z aitoff");

   // Draw the grid	
   TPad *pad2 = new TPad("pad2","",0,0,1,1);
   pad2->SetFillStyle(4000);
   pad2->SetFillColor(0);
   pad2->SetBorderSize(0);
   pad2->SetFrameBorderMode(0);
   pad2->SetFrameLineColor(0); 
   pad2->SetFrameBorderMode(0);
   pad2->Draw();
   pad2->cd();
   Double_t ymin = -89.5;
   Double_t ymax = 89.5;
   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

   pad2->Range(xmin-0.1*dx,ymin-0.1*dy,xmax+0.1*dx,ymax+0.1*dy);

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