#include "../Headers/draw.h"

int Backprojection(const Setup& config, const std::vector<Cone>& cones, 
                   std::vector<std::vector<double>>& image){
    // project cones onto the spherical surface
    std::cout << "Projecting cones onto the designated spherical surface..."<<std::endl;
    double progress=0.0;
    u_int currentConeNum(0);
    u_int totalConeNum(cones.size());
    double alpha(0);
    double beta(0);
    for (const Cone& event : cones)
    {
        // update progress bar
        int barWidth = 70;
        currentConeNum++;
        if (currentConeNum % 10 == 0)
        {
            progress = double(currentConeNum) / double(totalConeNum);
            std::cout << "[";
            int pos = barWidth * progress;
            for (int i = 0; i < barWidth; ++i) {
                if (i < pos) std::cout << "=";
                else if (i == pos) std::cout << ">";
                else std::cout << " ";
            }
            std::cout << "] " << int(progress * 100.0) << " %\r";
            std::cout.flush();
        }

        // update the image
        //TODO: add the uncertainty
        alpha = event.cosHalfAngle;
        for (int i = 0; i < config.thetaBins; i++)
        {
            for (int j = 0; j < config.phiBins; j++)
            {
                beta = getCosAngle(event.apex - config.xbs[i][j], event.axis);
                image[i][j] += std::exp(-(beta-alpha)*(beta-alpha) / 0.002);
            }
        }
    }
    return 0;
}

int drawImage(const Setup& config, const std::vector<std::vector<double>>& image){
    TH2D *histo = new TH2D("ROI", " ; Azimuth; Elevation", 
                          config.phiBins, -180, 180, 
                          config.thetaBins, -89.5, 89.5);
    TCanvas *canvas = new TCanvas("THCanvas","Sinocanvas", 1000, 500);
    for (int i = 1; i <= config.phiBins; i++)
    {
        for (int j = 1; j < config.thetaBins; j++)
        {
            histo->SetBinContent(i, j, image[j][i]);
        }
    }
    gStyle->SetOptStat(0);
    histo->GetZaxis()->SetLabelSize(0.02);
    histo->Draw("z aitoff");

    // grid
    float conv=M_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);
        }
    }
    // 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");
    // // // canvas1->cd(4)->SetLogz();

    // // // draw source 1
    // // TEllipse* el1 = new TEllipse(86.5276, 0, 5., 5.);
    // // el1->SetFillColor(0);
    // // el1->SetFillStyle(0);
    // // el1->SetLineColor(4);
    // // el1->Draw("SAME");

    // draw source 2
    // double sigma = 5;
    // TEllipse* el2 = new TEllipse(config.truePhi / conv, config.trueTheta / conv, sigma, sigma);
    // el2->SetFillColor(0);
    // el2->SetFillStyle(0);
    // el2->SetLineColor(4);
    // el2->Draw("SAME");

    // draw source in aitoff 
    double Ax, Ay;
    aitoff2xy(config.truePhi / conv, config.trueTheta / conv, Ax, Ay);
    TMarker *marker = new TMarker(Ax, Ay, 24);
    // marker->SetMarkerColorAlpha(kBlue, 0.0);
    marker->Draw("SAME");

    canvas->Draw();

    return 0;
}

int aitoff2xy(const double& l, const double& b, double &Al, double &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));
   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;

   return 0;
}