//-----------------------------------------------------------------
// Simplified version of plan view of Hall D beam dump.
// Script based on David's DrawGlueXBoundaries script
//
//-----------------------------------------------------------------

void Draw_HallD_BeamDump()
{
    // TColor::CreateColorWheel();

      char string[132];


      Double_t zmin = -600;
      Double_t zmax=2600;
      Double_t xmin= -500;
      Double_t xmax=1500;
      Double_t fudge=0.95;
      Double_t scale = fudge*(zmax-zmin)/(xmax-xmin);
       printf ("Scale factor =%d\n",scale);
       Double_t fe_density = 7.8;    // take Fe as cr shield
       Double_t dirt_density = 1.7;   // take dirt for overburden
       Double_t concrete_density = 2.7;   // concrete bunker

       Double_t xgrade=762;   // 25' grade above beamline
 
       Double_t dz=520;
       Double_t dr=dz/2;
       Double_t dc=91;
       Double_t dmuon=800;
       Double_t xmuon=136;

       Double_t zbdx=1600;
       Double_t Lbdx=250;
       // Double_t xbdx=110./2;     // note that in this confusing historical notation x is vertical, y is horizontal.
       // Double_t ybdx=205./2;    // note that in this confusing historical notation x is vertical, y is horizontal.
       Double_t xbdx=70./2;     // note that in this confusing historical notation x is vertical, y is horizontal.
       Double_t ybdx=165./2;    // note that in this confusing historical notation x is vertical, y is horizontal.

       Double_t dwall=30;
       Double_t dspace=30;

       Double_t zdump=150;
       Double_t Ldump=150;
       Double_t xdump=30;

       Double_t overburden_mwe=10*100;     // mwe converted to cm
       Double_t tshield= overburden_mwe/fe_density;
       Double_t zcrshield=zbdx-dspace;
       Double_t Lcrshield=Lbdx+2*dspace;
       Double_t xcrshield1=xgrade;
       Double_t xcrshield2=xgrade-tshield;

       Double_t xover2 = overburden_mwe/dirt_density + xbdx;
       Double_t arg = (xgrade-xbdx-dspace)/xover2 < 1 ? (xgrade-xbdx-dspace)/xover2 : 1;
       Double_t angle_over = acos(arg);
       Double_t zover = arg < 1 ? (zbdx+Lbdx/2) - sin(angle_over)*xover2 : 0;
       printf ("Assume: dirt density=%f, overburden=%f cm, angle=%f, zover = %f\n",dirt_density,xover2,angle_over*180/3.14159,zover);
  
      TCanvas *c1 = new TCanvas("c1","Draw_HallD_BeamDump",200,10,500*scale,500);

      TH2D *boundaries = new TH2D("boundaries", "Hall D Beam Dump Dimensions",zmax-zmin, zmin, zmax, xmax-xmin, xmin, xmax);
       boundaries->SetStats(0);
       boundaries->SetXTitle("Z (cm)");
       boundaries->SetYTitle("Y (cm)");
       boundaries->Draw();

       TLine *beamline = new TLine (zmin,0,zmax,0);
       beamline->SetLineStyle(4);
       beamline->Draw();

       TLine *grade = new TLine (zmin,xgrade,zmax,xgrade);
       grade->SetLineStyle(1);
       grade->Draw();

      sprintf (string,"Overburden = %.1f mwe\n",overburden_mwe/100);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.85,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      sprintf (string,"Dirt density = %.1f g/cm2\n",dirt_density);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.81,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      sprintf (string,"BDX length = %.1f cm\n",Lbdx);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.77,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      sprintf (string,"BDX height = 2x%.1f cm\n",xbdx);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.73,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      Double_t fe_shield = dz - (zdump + Ldump) +  dmuon; 
      Double_t concrete_shield = dc; 
      Double_t muon_shield = fe_shield + dc*concrete_density/fe_density; 

      sprintf (string,"Muon shield = %.1f cm Fe + %.1f cm Concrete\n",fe_shield,concrete_shield);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.69,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      sprintf (string,"Muon shield = %.1f cm Fe (effective)\n",muon_shield);
      printf("string=%s",string);
      t1 = new TLatex(0.12,0.65,string);
      t1->SetNDC();
      t1->SetTextSize(0.03);
      t1->Draw();

      // concrete surrounding

      const int Npoints=10;
      float xc[]={0,dr,dr,dr+dc,dr+dc,-dr-dc,-dr-dc,-dr,-dr,0.};
      float zc[]={dz,dz,-dz,-dz,dz+dc,dz+dc,-dz,-dz,dz,dz};

      TPolyLine *concrete = new TPolyLine(Npoints, zc, xc);
      concrete->SetLineWidth(2.0);
      concrete->SetLineColor(kBlack);
      concrete->Draw("l");
      concrete->SetFillStyle(3012);
      concrete->SetFillColor(kBlack);
      concrete->Draw("f");

      // muon iron

      const int Npoints=6;
      float xc[]={0,xmuon,xmuon,-xmuon,-xmuon,0};
      float zc[]={dz+dc,dz+dc,dz+dc+dmuon,dz+dc+dmuon,dz+dc,dz+dc};

      TPolyLine *dsteel = new TPolyLine(Npoints, zc, xc);
      dsteel->SetLineWidth(2.0);
      dsteel->SetLineColor(kBlack);
      dsteel->Draw("l");
      dsteel->SetFillStyle(3545);
      dsteel->SetFillColor(kBlack);
      dsteel->Draw("f");

      // steel casing

      const int Npoints=10;
      float xc[]={0,xdump,xdump,dr,dr,-dr,-dr,-xdump,-xdump,0.};
      float zc[]={zdump+Ldump,zdump+Ldump,0,0,dz,dz,0,0,zdump+Ldump,zdump+Ldump};

      TPolyLine *casing = new TPolyLine(Npoints, zc, xc);
      casing->SetLineWidth(2.0);
      casing->SetLineColor(kBlack);
      casing->Draw("l");
      casing->SetFillStyle(3545);
      casing->SetFillColor(kBlack);
      casing->Draw("f");

      // dump proper

      const int Npoints=6;
      float xc[]={0,xdump,xdump,-xdump,-xdump,0};
      float zc[]={zdump,zdump,zdump+Ldump,zdump+Ldump,zdump,zdump};

      TPolyLine *dump = new TPolyLine(Npoints, zc, xc);
      dump->SetLineWidth(2.0);
      dump->SetLineColor(kBlack);
      dump->Draw();

      TBox *filldump = new TBox(zdump,-xdump,zdump+Ldump,xdump);
      filldump->SetFillColor(590);
      filldump->Draw();

      // BDX experiment

      const int Npoints=6;
      float xc[]={0,xbdx,xbdx,-xbdx,-xbdx,0};
      float zc[]={zbdx,zbdx,zbdx+Lbdx,zbdx+Lbdx,zbdx,zbdx};

      TPolyLine *bdx = new TPolyLine(Npoints, zc, xc);
      bdx->SetLineWidth(2.0);
      bdx->SetLineColor(kBlack);
      bdx->Draw("l");
      bdx->SetFillColor(390);
      bdx->Draw("f");

      // BDX walls

      const int Npoints=9;
      float xc[]={0,-xbdx-dspace-dwall,-xbdx-dspace-dwall,0,0,-xbdx-dspace,-xbdx-dspace,0,0};
      float zc[]={zbdx-dspace-dwall,zbdx-dspace-dwall,zbdx+Lbdx+dspace+dwall,zbdx+Lbdx+dspace+dwall,zbdx+Lbdx+dspace,zbdx+Lbdx+dspace,zbdx-dspace,zbdx-dspace,zbdx-dspace-dwall};

      TPolyLine *bdxwall = new TPolyLine(Npoints, zc, xc);
      bdxwall->SetLineWidth(2.0);
      bdxwall->SetLineColor(kBlack);
      bdxwall->Draw("l");
      bdxwall->SetFillStyle(3012);
      bdxwall->SetFillColor(kBlack);
      bdxwall->Draw("f");

      for (int j=0; j<Npoints; j++) {
             xc[j] = -xc[j];
             }

      TPolyLine *bdxwall = new TPolyLine(Npoints, zc, xc);
      bdxwall->SetLineWidth(2.0);
      bdxwall->SetLineColor(kBlack);
      bdxwall->Draw("l");
      bdxwall->SetFillStyle(3012);
      bdxwall->SetFillColor(kBlack);
      bdxwall->Draw("f");

      // cosmic-ray shield experiment. Ignore if sufficient overburden already

     if (zover > 0) {

      const int Npoints=4;
      float xc[]={xgrade,xover2,xover2,xgrade,xgrade};
      float zc[]={zover,zbdx-dspace,zbdx+Lbdx+dspace,zbdx+Lbdx/2 + (zbdx+Lbdx/2-zover),zover};

      TPolyLine *crover = new TPolyLine(Npoints, zc, xc);
      crover->SetLineWidth(2.0);
      crover->SetLineColor(kBlack);
      crover->Draw("l");
      crover->SetFillColor(kBlack);
      crover->SetFillStyle(3003);
      crover->Draw("f");
      }


      sprintf (string,"Draw_HallD_BeamDump.pdf");
      c1->SaveAs(string);
      sprintf (string,"Draw_HallD_BeamDump.png");
      c1->SaveAs(string);

}