#include <iostream>
#include <string>
#include <TROOT.h>
#include <TChain.h>
#include <TTree.h>
#include <TH1D.h>

void doseplots(int id=0, const char *var=0, const char* sel=0, TH1 *hproj=0)
{
   gSystem->Load("deq_ambient_C.so");
   gSystem->Load("h1dose_C.so");

   TChain *h1 = new TChain("h1");
   TChain *h2 = new TChain("h2");
   TChain *h3 = new TChain("h3");
   TChain *h4 = new TChain("h4");
   TChain *h5 = new TChain("h5");
   TChain *h6 = new TChain("h6");

   for (int seg=8; seg < 9; ++seg) {
      char rootf[50];
      sprintf(rootf, "results_%d.root", seg*10000);
      h1->Add(rootf);
      h2->Add(rootf);
      h3->Add(rootf);
      h4->Add(rootf);
      h5->Add(rootf);
      h6->Add(rootf);
   }

   int microcounts = h1->Draw("pin[4]","kind==3&&pin[4]>2.8&&desum>2e-3&&pin[4]>3.0&&pin[4]<3.6");
   double duration = microcounts / 250e6;
   std::cout << "simulation covers " << duration 
             << " s at nominal intensity." << std::endl;

   // dose plot inside shielded tagger readout box
   if (h4 == 0) {
      std::cout << "h4 not found!" << std::endl;
      return;
   }
   h1dose *s4 = new h1dose(h4,1,-212.,-132.,0,-75.,-60.,2,710.,860.,duration);
   h4->Process(s4);
   TH1D *d4 = (TH1D*)gROOT->FindObject("h4_dose");
   if (d4) {
      d4->SetMinimum(0);
      d4->SetTitle("neutrons inside TAGM readout shielding box");
      d4->GetXaxis()->SetTitle("z (cm)");
      d4->GetYaxis()->SetTitle("dose (mrem/hr)");
      d4->Draw("E");
   }
   else
      std::cout << "h4_dose not found!" << std::endl;

   // dose plot along y at tagger microscope on focal plane
   if (h5 == 0) {
      std::cout << "h5 not found!" << std::endl;
      return;
   }
   h1dose *s5y = new h1dose(h5,0,-500.,0.,10,-180.,180.,1,710.,860.,duration);
   h5->Process(s5y);
   TH1D *d5y = (TH1D*)gROOT->FindObject("h5_dose");
   d5y->SetName("h5_dosey");
   if (d5y) {
      d5y->SetMinimum(0);
      d5y->SetTitle("neutrons at focal plane averaged along z");
      d5y->GetXaxis()->SetTitle("y from tagger midplane (cm)");
      d5y->GetYaxis()->SetTitle("dose (mrem/hr)");
      d5y->Draw("E");
   }
   else
      std::cout << "h5_dosey not found!" << std::endl;

   // dose plot along tagger focal plane
   if (h5 == 0) {
      std::cout << "h5 not found!" << std::endl;
      return;
   }
   h1dose *s5 = new h1dose(h5,0,-500.,0.,1,-180.,180.,10,415.,1258.,duration);
   h5->Process(s5);
   TH1D *d5 = (TH1D*)gROOT->FindObject("h5_dose");
   if (d5) {
      d5->SetMinimum(0);
      d5->SetTitle("neutrons on focal plane averaged over y");
      d5->GetXaxis()->SetTitle("z (cm)");
      d5->GetYaxis()->SetTitle("dose (mrem/hr)");
      d5->Draw("E");
   }
   else
      std::cout << "h5_dose not found!" << std::endl;

   // dose plot near the floor under the tagger
   if (h6 == 0) {
      std::cout << "h6 not found!" << std::endl;
      return;
   }
   h1dose *s6 = new h1dose(h6,1,-275.,275.,0,-500.,500.,10,325.,875.,duration);
   h6->Process(s6);
   TH1D *d6 = (TH1D*)gROOT->FindObject("h6_dose");
   if (d5) {
      d6->SetMinimum(0);
      d6->SetTitle("neutrons near floor under tagger magnet");
      d6->GetXaxis()->SetTitle("z (cm)");
      d6->GetYaxis()->SetTitle("dose (mrem/hr)");
      d6->Draw("E");
   }
   else
      std::cout << "h6_dose not found!" << std::endl;

   if (hproj == 0) {
      if (id == 1 && h1)
         h1->Draw(var,sel);
      else if (id == 2 && h2)
         h2->Draw(var,sel);
      else if (id == 3 && h3)
         h3->Draw(var,sel);
      else if (id == 4 && h4)
         h4->Draw(var,sel);
      else if (id == 5 && h5)
         h5->Draw(var,sel);
      else if (id == 6 && h6)
         h6->Draw(var,sel);
   }
   else {
      if (id == 1 && h1)
         h1->Project(hproj->GetName(),var,sel);
      else if (id == 2 && h2)
         h2->Project(hproj->GetName(),var,sel);
      else if (id == 3 && h3)
         h3->Project(hproj->GetName(),var,sel);
      else if (id == 4 && h4)
         h4->Project(hproj->GetName(),var,sel);
      else if (id == 5 && h5)
         h5->Project(hproj->GetName(),var,sel);
      else if (id == 6 && h6)
         h6->Project(hproj->GetName(),var,sel);
   }
}

void Espectrum(int id, TH1D *hspec)
{
   std::string title;
   title = std::string(hspec->GetTitle()) + " photon hits";
   TH1D *hphot = hspec->Clone("hphot");
   hphot->SetTitle(title.c_str());
   doseplots(id,"pin[3]","kind==1",hphot);
   hphot->SetFillColor(5);

   title = std::string(hspec->GetTitle()) + " electron hits";
   TH1D *helec = hspec->Clone("helec");
   helec->SetTitle(title.c_str());
   doseplots(id,"pin[3]","kind==3",helec);
   helec->SetFillColor(2);

   title = std::string(hspec->GetTitle()) + " neutron hits";
   TH1D *hneut = hspec->Clone("hneut");
   hneut->SetTitle(title.c_str());
   doseplots(id,"pin[3]","kind==13",hneut);
   hneut->SetFillColor(7);
   
   title = std::string(hspec->GetTitle()) + " hits";
   THStack *hs = new THStack("hs",title.c_str());
   hs->Add(hphot);
   hs->Add(helec);
   hs->Add(hneut);
   hs->Draw();
   hs->GetHistogram()->GetXaxis()->SetTitle("particle energy (GeV)");
   hs->GetHistogram()->GetYaxis()->SetTitle("counts");
   hs->Draw();
}

void spectrum(int id, TH1D *hspec, const char* var)
{
   std::string title;
   title = std::string(hspec->GetTitle()) + " photon hits";
   TH1D *hphot = hspec->Clone("hphot");
   hphot->SetTitle(title.c_str());
   doseplots(id,var,"kind==1",hphot);
   hphot->SetFillColor(5);

   title = std::string(hspec->GetTitle()) + " electron hits";
   TH1D *helec = hspec->Clone("helec");
   helec->SetTitle(title.c_str());
   doseplots(id,var,"kind==3",helec);
   helec->SetFillColor(2);

   title = std::string(hspec->GetTitle()) + " neutron hits";
   TH1D *hneut = hspec->Clone("hneut");
   hneut->SetTitle(title.c_str());
   doseplots(id,var,"kind==13",hneut);
   hneut->SetFillColor(7);
   
   title = std::string(hspec->GetTitle()) + " hits";
   THStack *hs = new THStack("hs",title.c_str());
   hs->Add(hphot);
   hs->Add(helec);
   hs->Add(hneut);
   hs->Draw();
   hs->GetHistogram()->GetXaxis()->SetTitle("energy deposition (GeV)");
   hs->GetHistogram()->GetYaxis()->SetTitle("counts");
   hs->Draw();
}