#include <TROOT.h>
#include <TH2D.h>
#include "MetroProMap.h"
#include "Map2D.h"
#include "ablator.h"

const double thin_dx = 5.8;
const double thin_dy = 5.8;

const double diam_dx = 6.0;
const double diam_dy = 6.0;

double beamspot(double sigma_x=1.200, double sigma_y=0.800)
{

   // form Map2D images of the beam and diamond target
   Map2D *view = new Map2D(TH2D("view","viewport",1000,-5,5,1000,-5,5));
   Map2D *diam = new Map2D(*view);
   Map2D *beam = new Map2D(*view);
   diam->SetName("diam");
   beam->SetName("beam");
   for (int ix=1; ix<=1000; ++ix) {
      for (int iy=1; iy<=1000; ++iy) {
         double x = view->GetXaxis()->GetBinCenter(ix);
         double y = view->GetYaxis()->GetBinCenter(iy);
         double bintens = exp(-0.5*(pow(x/sigma_x,2)+pow(y/sigma_y,2)));
         double dprofil = (fabs(x+y)/sqrt(2) < thin_dx/2 &&
                           fabs(x-y)/sqrt(2) < thin_dy/2)?  20 : 500;
         beam->SetBinContent(ix,iy,bintens);
         diam->SetBinContent(ix,iy,dprofil);
      }
   }

   // smooth the ablated edges of the diamond target window frame
   ablator abl;
   abl.smooth(*diam);
   Map2D *diam_smoothed = (Map2D*)gROOT->FindObject("diam_smoothed");
   *diam = *diam_smoothed;
   diam->SetName("diam");
   for (int ix=1; ix<=1000; ++ix) {
      for (int iy=1; iy<=1000; ++iy) {
         double x = view->GetXaxis()->GetBinCenter(ix);
         double y = view->GetYaxis()->GetBinCenter(iy);
         double dprofil = (fabs(x+y)/sqrt(2) < diam_dx/2 &&
                           fabs(x-y)/sqrt(2) < diam_dy/2)?
                          diam->GetBinContent(ix,iy) : 0;
         diam->SetBinContent(ix,iy,dprofil);
      }
   }

   // find the beam intensity averaged target thickness
   Map2D *prod = new Map2D(*view);
   prod->SetName("prod");
   prod->SetTitle("beam * target profile product");
   diam->SetTitle("diamond target profile");
   beam->SetTitle("beam intensity profile");
   prod->Multiply(beam,diam);
   return prod->Integral()/beam->Integral();
}