#include <iostream>
#include <math.h>

#include <TROOT.h>
#include <TH2D.h>
#include <TCanvas.h>

TCanvas *c1=0;

TH2D *mounts()
{
   // Geometry description:
   // The mounted diamond is described by the following geometry layout:
   //
   //         +-------------------------------+
   //         |  ^                  Al(t4)    |
   //         |  |           w4               |
   //         |<-+--------------------------->|
   //         |  |   +-----------------+ ^    |
   //         |  |   |                 | |    |
   //         |  |   |   ###########   | |    |
   //         |  |   |   # diamond #   | |    |
   //         |  |   |\  #   (t1)  #   | | h3 |
   //         |  |   |\\\#         #   | |    |
   //         |  |   |\\\#####(t2)##   | |    |
   //         |  |   |\\(t3)\          | |    |
   //         |  |   +-----------------+ v    |
   //         |  |   <----------------->      |
   //         |  | h4        w3               |
   //         |  v                            |
   //         +-------------------------------+
   //
   // The diamond itself is shown by the hashed outline in the figure.
   // It has an inner target area of thickness t1 and dimensions w1 x h1.
   // This area is centered inside the overall diamond outline of size
   // w2 x h2, and outer frame thickess t2. The diamond is centered inside
   // a rectangular opening in the aluminum mount of dimensions w3 x h3,
   // which is completely open except for a thin triangular membrane that
   // covers the lower left corner of the opening and holds the diamond
   // target by its corner. The thickness of the membrane is t3. The
   // aluminum mount around the target window is of thickness t4 and has
   // outer dimensions no less than w4 x h4. All dimensions are given in
   // meters.
   //
   // The output from the calculation is a map of the above geometry where
   // every pixel represents the relative intensity of the bremsstrahlung
   // generated by the electron beam at that point. The electron beam is
   // represented by a Gaussian plus a uniform halo.

   double w1=6.8e-3, h1=6.8e-3, t1=20e-6;
   double w2=7.3e-3, h2=7.3e-3, t2=300e-6;
   double w3=12.0e-3, h3=12.0e-3, t3=800e-6;
   double w4=15.0e-3, h4=15.0e-3, t4=3e-3;

   double radlen_diamond=15e-2;
   double radlen_aluminum=8.9e-2;

   double beam_sigmaX=1.1e-3;
   double beam_sigmaY=0.9e-3;
   double beam_halo_fraction=1e-5;

   double image_width = w4;
   double image_height = h4;

   int pixdim = 500;
   TH2D *rintens = new TH2D("rintens","radiated intensity",
                            pixdim,-image_width/2,+image_width/2,
                            pixdim,-image_height/2,+image_height/2);

   double beam_halo_norm = 0;
   double beam_cent_norm = 0;
   for (int ix=1; ix <= pixdim; ++ix) {
      double x = rintens->GetXaxis()->GetBinCenter(ix);
      for (int iy=1; iy <= pixdim; ++iy) {
         double y = rintens->GetYaxis()->GetBinCenter(iy);
         double bd[2] = {x/beam_sigmaX, y/beam_sigmaY};
         beam_cent_norm += exp(-0.5*(bd[0]*bd[0]+bd[1]*bd[1]));
         beam_halo_norm += 1;
      }
   }

   double active_sum=0;
   double frame_sum=0;
   double alum_sum=0;
   for (int ix=1; ix <= pixdim; ++ix) {
      double x = rintens->GetXaxis()->GetBinCenter(ix);
      for (int iy=1; iy <= pixdim; ++iy) {
         double y = rintens->GetYaxis()->GetBinCenter(iy);
         double bd[2] = {x/beam_sigmaX, y/beam_sigmaY};
         double bintens = beam_halo_fraction/beam_halo_norm
                        + exp(-0.5*(bd[0]*bd[0]+bd[1]*bd[1]))/
                          beam_cent_norm;
         double radthick;
         if (fabs(x) < w1/2 && fabs(y) < h1/2) {
            radthick = t1/radlen_diamond;
            active_sum += radthick*bintens;
         }
         else if (fabs(x) <  w2/2 && fabs(y) < h2/2) {
            radthick = t2/radlen_diamond;
            frame_sum += radthick*bintens;
         }
         else if (fabs(x) >  w3/2 || fabs(y) > h3/2) {
            radthick = t4/radlen_aluminum;
            alum_sum += radthick*bintens;
         }
         else if (x+y < 1e-3-(w1+w2)/2) {
            radthick = t3/radlen_aluminum;
            alum_sum += radthick*bintens;
         }
         else {
            radthick = 0;
         }
         rintens->SetBinContent(ix,iy,radthick*bintens);
      }
   }
   std::cout << "active, frame, aluminum yields are: "
             << active_sum << ", "
             << frame_sum << ", "
             << alum_sum << std::endl;
   double total_sum = active_sum+frame_sum+alum_sum;
   std::cout << "active, frame, aluminum factions are: "
             << active_sum/total_sum << ", "
             << frame_sum/total_sum << ", "
             << alum_sum/total_sum << std::endl;

   c1 = (TCanvas*)gROOT->FindObject("c1");
   if (c1) {
      c1->cd();
   }
   else {
      c1 = new TCanvas("c1","c1",10,10,550,500);
   }
   c1->SetLeftMargin(0.10);
   c1->SetRightMargin(0.22);
   c1->SetTopMargin(0.10);
   rintens->SetContour(100);
   rintens->SetStats(0);
   rintens->Draw("colz");
   c1->Update();
   return rintens;
}