////////////////////////////////////////////////////
///
/// Plot Target Related Quantities for Beam Monitoring
///
/// Authors: Tom Osiecki + Brett Biren
///
/// osiecki@mail.hep.utexas.edu
/// bv@bnl.gov
///
////////////////////////////////////////////////////

#include "TargetModule.h"

#include "StripHist.h"

#include <JobControl/JobCModuleRegistry.h>
#include <JobControl/JobCResult.h>
#include <MessageService/MsgService.h>

#include <RawData/RawRecord.h>
#include <RawData/RawBeamMonBlock.h>
#include <RawData/RawBeamData.h>

#include <Conventions/Munits.h>

#include <BeamDataUtil/BDTarget.h>
#include <BeamDataUtil/BDSwicCalibrator.h>

#include <TGraph.h>
#include <TGaxis.h>
#include <TCanvas.h>
#include <TLatex.h>
#include <TEllipse.h>
#include <TBox.h>

//ClassImp(TargetModule)

#include <cmath>

#define minimum(a,b) ((a) > (b) ? (b) : (a)) 
#define maximum(a,b) ((a) > (b) ? (a) : (b)) 

using namespace std;

CVSID("$Id: TargetModule.cxx,v 1.14 2006/05/27 07:28:39 rhatcher Exp $");
JOBMODULE(TargetModule,"MonTgt","Target quantities related for Monitoring");

TargetModule::TargetModule()
{
  StripHist *sh1 = new StripHist("Beam Size H","Beam Size Red=Horizontal Blue=Vertical");
  StripHist *sh5 = new StripHist("Beam Size V","Beam Size Red=Horizontal Blue=Vertical");
  StripHist *sh2 = new StripHist("Fraction Target","The % of the beam that hit the target");
  StripHist *sh3 = new StripHist("Fraction Baffle","The % of the beam that hit the baffle");
  StripHist *sh4 = new StripHist("Integral Flux NuMI","Integral flux of the beam in 10^{12} units Red=NuMI Blue=Target");
  StripHist *sh6 = new StripHist("Integral Flux Tgt","Integral flux of the beam in 10^{12} units Red=NuMI Blue=Target");
  
  TH2F *BeamParameters = new TH2F("Beam Parameters","Red Contours Represent 1,2,3 #sigma of Beam",10,-10,10,10,-10,10);

  BeamParameters->GetXaxis()->SetTitle("X (mm)");
  BeamParameters->GetYaxis()->SetTitle("Y (mm)");
  BeamParameters->GetXaxis()->CenterTitle();
  BeamParameters->GetYaxis()->CenterTitle();

  sh1->SetStripRange(1*Munits::day);
  sh2->SetStripRange(1*Munits::day);
  sh3->SetStripRange(1*Munits::day);
  sh4->SetStripRange(1*Munits::day);
  sh5->SetStripRange(1*Munits::day);
  sh6->SetStripRange(1*Munits::day);

  sh1->SetLineColor(2);
  sh5->SetLineColor(4);
  sh4->SetLineColor(2);
  sh6->SetLineColor(4);

  sh1->SetGraphYTitle("Beam Size (mm)");
  sh5->SetGraphYTitle("Beam Size (mm)");
  sh2->SetGraphYTitle("Fraction on Target (%)");
  sh3->SetGraphYTitle("Fraction on Baffle (%)");
  sh4->SetGraphYTitle("Integral Flux (10^{12} protons)");
  sh6->SetGraphYTitle("Integral Flux (10^{12} protons)");
  
  fStripHist["Beam Size H"]        = sh1;
  fStripHist["Beam Size V"]        = sh5;
  fStripHist["Fraction Target"]    = sh2;
  fStripHist["Fraction Baffle"]    = sh3;
  fStripHist["Integral Flux NuMI"] = sh4; 
  fStripHist["Integral Flux Tgt"]  = sh6;

  fTH2F["Beam Parameters"] = BeamParameters;
}

TargetModule::~TargetModule()
{

}

void TargetModule::BeginJob()
{
  HistMan hm = this->GetHistMan();

  StripHist *sh1 = fStripHist["Beam Size H"];
  StripHist *sh5 = fStripHist["Beam Size V"];
  
  TCanvas* canvas1 = new TCanvas("Beam Size","Beam Size",500,400);

  sh1->DrawStrip("AL");
  sh5->DrawStrip("LP");

  hm.Adopt("Target",canvas1);
  
  StripHist *sh2 = fStripHist["Fraction Target"];
  
  TCanvas *canvas2 = new TCanvas("Fraction Beam on Target",sh2->GetHist().GetTitle(),500,400);

  sh2->DrawStrip("AL");
  
  hm.Adopt("Target",canvas2);

  StripHist *sh3 = fStripHist["Fraction Baffle"];
  
  TCanvas *canvas3 = new TCanvas("Fraction Beam on Baffle",sh3->GetHist().GetTitle(), 500,400);

  sh3->DrawStrip("AL");

  hm.Adopt("Target",canvas3);

  StripHist *sh4 = fStripHist["Integral Flux NuMI"];
  StripHist *sh6 = fStripHist["Integral Flux Tgt"];

  TCanvas *canvas4 = new TCanvas("Integral Flux","Integral Flux",500,400);

  sh4->DrawStrip("AL");
  sh6->DrawStrip("LP");

  hm.Adopt("Target",canvas4);

  //TH2F *th1 = fTH2F["Beam Parameters"];

  TCanvas *canvas5 = new TCanvas("Target Summary Canvas","Target Summary Canvas",600,600);
  canvas5->Divide(2,3);
  canvas5->cd(1);
  //Done in Fill();
  canvas5->cd(2);
  sh1->DrawStrip("AL");
  sh5->DrawStrip("LP");
  canvas5->cd(3);
  //Done in Fill();
  canvas5->cd(4);
  sh2->DrawStrip("AL");
  canvas5->cd(5);
  sh4->DrawStrip("AL");
  sh6->DrawStrip("LP");
  canvas5->cd(6);
  sh3->DrawStrip("AL");

  hm.Adopt("Summary Canvases",canvas5);

}

static double approx_gaussian(double x, double sigma)
{
    //Returns an approximate integral of a gaussian centered at (0,0) from [0, x].
    const double sqrt2 = sqrt(2.0);
    
    return 0.5 - 0.5 * ((sigma * sqrt2) / (3.14159 * 2.0)) * exp(-(x) / (sqrt2 * sigma));
}

static double percentageOnTarget(double x, double y, double sigmax, double sigmay)
{
  if( (pow(x+1.5,2) + pow(y-1.0,2)) < 5.5 && x > -4.7 && x < 1.7) {
    double answerx = approx_gaussian(1.7, sigmax) + approx_gaussian(4.7, sigmax);
    double miny    = sqrt(5.5-pow(x+1.5,2))+1.0;
    double maxy    = sqrt(5.5-pow(x+1.5,2))+1.0;
    double answery = approx_gaussian(miny, sigmay) + approx_gaussian(maxy, sigmay);
    return answerx * answery;
  } else {
    return 0.0;
  }
}

static double percentageOnBaffle(double percentage_on_target)
{
  if(percentage_on_target < 1.0) {
    return 1.0 - percentage_on_target;
  } else {
    return 0.0;
  }
}

void TargetModule::Fill(const RawBeamMonHeaderBlock& head, const RawBeamMonBlock& block)
{
  MSG("BD",Msg::kDebug) << "Starting TargetModule::Fill" << endl;

  //BDSwicCalibrator::Get().Calibrate(head,block);

  const RawBeamData *Emtgths = block["E:MTGTHS"]; // Gives Horizontal Beam Sigma
  const RawBeamData *Emtgtvs = block["E:MTGTVS"]; // Gives Vertical Beam Sigam  
  const RawBeamData *Ehp121  = block["E:HP121"];  
  const RawBeamData *Ehptgt  = block["E:HPTGT"];
  const RawBeamData *Evp121  = block["E:VP121"];
  const RawBeamData *Evptgt  = block["E:VPTGT"];
  const RawBeamData *Etortgt = block["E:TORTGT"];
  const RawBeamData *Etrtgtd = block["E:TRTGTD"];
  //  const RawBeamData *Inutarz = block["I:NUTARZ"];

  //if(Inutarz) MSG("BD",Msg::kDebug) << "config=" << Inutarz->GetData()[0] << endl;

  if(!Emtgths || !Emtgtvs || !Ehp121 || !Ehptgt || !Evp121 || !Evptgt || !Etortgt || !Etrtgtd) {
    MSG("BD",Msg::kDebug) << "Do not have the devices we need!" << endl;
    if(!Emtgths) MSG("BD",Msg::kDebug) << "  No MTGTHS" << endl;
    if(!Emtgtvs) MSG("BD",Msg::kDebug) << "  No MTGTVS" << endl;
    if(!Ehp121)  MSG("BD",Msg::kDebug) << "  No HP121" << endl;
    if(!Ehptgt)  MSG("BD",Msg::kDebug) << "  No HPTGT" << endl;
    if(!Evp121)  MSG("BD",Msg::kDebug) << "  No VP121" << endl;
    if(!Evptgt)  MSG("BD",Msg::kDebug) << "  Np VPTGT" << endl;
    if(!Etortgt) MSG("BD",Msg::kDebug) << "  No TORTGT" << endl;
    if(!Etrtgtd) MSG("BD",Msg::kDebug) << "  No TRTGTD" << endl;
   return;
  }

  if(!(Emtgths -> GetDataLength() &&
       Emtgtvs -> GetDataLength() && 
       Ehp121  -> GetDataLength() &&
       Ehptgt  -> GetDataLength() &&
       Evp121  -> GetDataLength() &&
       Evptgt  -> GetDataLength() &&
       Etortgt -> GetDataLength() && 
       Etrtgtd -> GetDataLength()))
    return;
    

  StripHist* sh1 = fStripHist["Beam Size H"];
  StripHist* sh5 = fStripHist["Beam Size V"];
  StripHist* sh2 = fStripHist["Fraction Target"];
  StripHist* sh3 = fStripHist["Fraction Baffle"];
  StripHist* sh4 = fStripHist["Integral Flux NuMI"];
  StripHist *sh6 = fStripHist["Integral Flux Tgt"];

  TH2F *th1 = fTH2F["Beam Parameters"];

  double sigmah = Emtgths->GetData()[0];
  double sigmav = Emtgtvs->GetData()[0];
  double dae    = Etortgt->GetSeconds() + 1.0e-6*Etortgt->GetMsecs();
  double hp121  = Ehp121->GetData()[0];
  double hptgt  = Ehptgt->GetData()[0];
  double vp121  = Evp121->GetData()[0];
  double vptgt  = Evptgt->GetData()[0];
  double P2tgt  = Etrtgtd->GetData()[0];
  double pos    = 1168.04; //For LE.  ME=1164.8 HE=1159.8

  if(P2tgt < 0.1) {
    //If Beam intensity low, get out.
    return;
  }

  double X = hptgt + (((hptgt - hp121) / (1134.9 - 1095.8)) * (pos - 1134.9));
  double Y = vptgt + (((vptgt - vp121) / (1135.8 - 1096.7)) * (pos - 1135.8));
  
  MSG("BD",Msg::kDebug) << "Beam Pos : X=" << X << " +/- " << sigmav << " Y=" << Y << " +/- " << sigmah << endl;
  MSG("BD",Msg::kDebug) << "   Also, hptgt=" << hptgt << " : vptgt=" << vptgt << " : hp121=" << hp121 << " : vp121=" << vp121 << endl;
  

  /* For when this works */
  BDTarget *target = new BDTarget();
  target->SetSpill(head, block);
  double X2 = 0.0, Y2 = 0.0, sigmav2 = 0.0, sigmah2 = 0.0;
  int blah = 0;
  blah = target->ProfileProjection(X2,Y2,sigmav2,sigmah2);
  MSG("BD",Msg::kDebug) << "Beam Pos : X2=" << X2<< " +/- " << sigmav2 << " Y2=" << Y2 << " +/- " << sigmah2 
			<< " and blah=" << blah << endl;
  

  sh1->Fill(dae, sigmah);
  sh5->Fill(dae, sigmav);

  double min1 = 0.0, max1 = 0.0;

  RangeFinder(sh1, sh5, min1, max1);
  
  sh1->GetStrip().GetYaxis()->SetRangeUser(min1, max1);
  sh5->GetStrip().GetYaxis()->SetRangeUser(min1, max1);

  double PerOnTarget = percentageOnTarget(X,Y,sigmav,sigmah);
  double PerOnBaffle = percentageOnBaffle(PerOnTarget);

  sh2->Fill(dae, PerOnTarget*100.0);
  sh3->Fill(dae, PerOnBaffle*100.0);

  double min2 = 0.0, max2 = 0.0;

  RangeFinder(sh2, sh3, min2, max2);

  sh2->GetStrip().GetYaxis()->SetRangeUser(min2, max2);
  sh3->GetStrip().GetYaxis()->SetRangeUser(min2, max2);

  HistMan hm = this->GetHistMan();
  TCanvas *can = (TCanvas*) hm.GetObject("Summary Canvases/Target Summary Canvas");
  TVirtualPad *pad = can->cd(3);
  pad->Clear();

  TLatex l;
  l.SetTextSize(0.1);  
  l.DrawLatex(0.3,0.8,"Beam Parameters");
  l.DrawLatex(0.1,0.5,Form("X=%.4f mm",X));
  l.DrawLatex(0.5,0.5,Form("#sigma_{X}=%.4f mm",sigmav));
  l.DrawLatex(0.1,0.2,Form("Y=%.4f mm",Y));
  l.DrawLatex(0.5,0.2,Form("#sigma_{Y}=%.4f mm",sigmah));


  TVirtualPad *pad2 = can->cd(1);
  pad2->Clear();
  th1->SetStats(0);
  th1->SetTitle("");
  th1->Draw();
  TEllipse *el = new TEllipse(-1.0,1.0,5.5,5.5);
  el->SetLineWidth(2);
  el->SetLineColor(13);
  el->Draw();
  TBox *box = new TBox(-4.7, -6.5, 1.7, 8.5);
  box->SetFillStyle(0);
  box->SetLineWidth(2);
  box->SetLineColor(13);
  box->Draw();
  //Now Draw Beam
  TEllipse *e2 = new TEllipse(X,Y,sigmav,sigmah);
  e2->SetLineColor(2);
  e2->Draw();
  TEllipse *e3 = new TEllipse(X,Y,2*sigmav,2*sigmah);
  e3->SetLineColor(2);
  e3->Draw();
  TEllipse *e4 = new TEllipse(X,Y,3*sigmav,3*sigmah);
  e4->SetLineColor(2);
  e4->Draw();
  
  pad2->Update();

  int N = sh4->GetStrip().GetN();
  if(N == 0) {
    sh4->Fill(dae, P2tgt);
    sh6->Fill(dae, P2tgt * PerOnTarget);
  } else {
    double x1 = 0.0, y1 = 0.0;
    sh4->GetStrip().GetPoint(N-1, x1, y1);
    sh4->Fill(dae, (P2tgt + y1));
    double x2 = 0.0, y2 = 0.0;
    sh6->GetStrip().GetPoint(N-1, x2, y2);
    sh6->Fill(dae, (P2tgt * PerOnTarget + y2));
  }

  double min3 = 0.0, max3 = 0.0;
  
  RangeFinder(sh4, sh6, min3, max3);
  
  sh4->GetStrip().GetYaxis()->SetRangeUser(min3, max3);
  sh6->GetStrip().GetYaxis()->SetRangeUser(min3, max3);
  //sh4->GetStrip().GetYaxis()->SetNoExponent(1);
  //sh6->GetStrip().GetYaxis()->SetNoExponent(1);

}

void TargetModule::RangeFinder(StripHist *sh1, StripHist *sh2, double &min, double &max)
{
  double max1 = sh1->GetMax(), max2 = sh2->GetMax();
  double min1 = sh1->GetMin(), min2 = sh2->GetMin();

  min = minimum(min1, min2);
  max = maximum(max1, max2);

  if(fabs(min) > 0.001) {
    if(min < 0) {
      min *= 1.5;
    } else {
      min *= 0.5;
    }
  } else {
    min = -0.5;
  }

  if(fabs(max) > 0.001) {
    if(max < 0) {
      max *= 0.5;
    } else {
      max *= 1.5;
    }
  } else {
    max = 0.5;
  }


}