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

#include "LossModule.h"

#include "StripHist.h"

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

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

#include <DataUtil/GetRecords.h>
#include <BeamDataUtil/BDDevices.h>

#include <MessageService/MsgService.h>

#include <Conventions/Munits.h>

#include <TGraph.h>
#include <TCanvas.h>
#include <TH2F.h>
#include <TLatex.h>
#include <TStyle.h>
#include <TGaxis.h>

#include <vector>

//ClassImp(LossModule)

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

using namespace std;

CVSID("$Id: LossModule.cxx,v 1.9 2005/06/09 15:26:43 thosieck Exp $");
JOBMODULE(LossModule,"MonLoss","Proton loss related for Monitoring");

LossModule::LossModule()
{
 

  for (int bafn = 1; bafn <= 2; ++bafn) {
    
    StripHist* sh = new StripHist(Form("Baffle Scrape %d",bafn),
				  Form("Estimate of beam loss on baffle, thermocouple %d",bafn),
				  100,-3000,100);
    sh->SetStripRange(1*Munits::day);
    sh->GetHist().SetXTitle("(Baffle Temp - Air Temp)/Toroid (degC/1E12 protons)");
    fStripHist[Form("Baffle Scrape %d",bafn)] = sh;
  }

  TGraph *sh1 = new TGraph("BPM Location along Beam Line H=Red V=Blue","BPM Location along Beam Line H=Red V=Blue");
  fTGraph["Beam Position H"] = sh1;
  TGraph *sh2 = new TGraph("BPM Location along Beam Line H=Red V=Blue","BPM Location along Beam Line H=Red V=Blue");
  fTGraph["Beam Position V"] = sh2;
  
  sh1->GetXaxis()->SetTitle("Distance from Q806 (ft)");
  sh2->GetXaxis()->SetTitle("Distance from Q806 (ft)");

  TH1F *s1 = new TH1F("Beam Loss","Beam Loss along tunnel",100,-40,1200);

  s1->SetXTitle("Distance from Q806 (ft)");
  s1->SetYTitle("Beam Loss (Rad / s)");
  s1->GetXaxis()->CenterTitle();
  s1->GetYaxis()->CenterTitle();
  s1->GetYaxis()->SetTitleOffset(1.3);
  //TGaxis::SetMaxDigits(1);

  fTH1F["Beam Loss"] = s1;
}

LossModule::~LossModule()
{

}

typedef std::map<std::string,StripHist*>  StripHistMap;
typedef std::map<std::string,TH1F*> TH1FMap;

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

  StripHist *sht1 = fStripHist["Baffle Scrape 1"];
  TCanvas* canvas = new TCanvas("Baffle Scrape 1",sht1->GetHist().GetTitle(),500,400);
  sht1->Draw("AL");
  hm.Adopt("Loss",canvas);

  StripHist *sht2 = fStripHist["Baffle Scrape 2"];
  TCanvas* canvas2 = new TCanvas("Baffle Scrape 2",sht2->GetHist().GetTitle(),500,400);
  sht2->Draw("AL");
  hm.Adopt("Loss",canvas2);

  TGraph *sht3 = fTGraph["Beam Position H"];
  TGraph *sht4 = fTGraph["Beam Position V"];

  sht3->SetLineColor(2);
  sht4->SetLineColor(4);

  TCanvas *canvas3 = new TCanvas("Beam Position","BPM Location along Beam Line",500,400);
  sht3->Draw("AL");
  sht4->Draw("LP");
  hm.Adopt("Loss",canvas3);

  TH1FMap::iterator mit2, done2=fTH1F.end();
  for(mit2=fTH1F.begin(); mit2 != done2; ++mit2) {
    TH1F *sh2 = mit2->second;
    TCanvas *canvas3 = new TCanvas(sh2->GetTitle(),sh2->GetTitle(),500,400);
    TVirtualPad* pad = gPad, *oldpad = gPad;    pad->cd();
    sh2->Draw();
    oldpad->cd();
    hm.Adopt("Loss",canvas3);
  }

  TH1F *th1 = fTH1F["Beam Loss"];

  TCanvas *canvas4 = new TCanvas("Proton Beam Transport Line Summary Canvas",
				 "Proton Beam Transport Line Summary Canvas",800,800);
  canvas4->Divide(1,3);
  canvas4->cd(1);
  sht3->Draw("AL");
  sht4->Draw("LP");
  canvas4->cd(2);
  th1->Draw();
  canvas4->cd(3);
  
  hm.Adopt("Summary Canvases",canvas4);

}

void LossModule::Fill(const RawBeamMonHeaderBlock& /*head*/, const RawBeamMonBlock& block)
{
  const RawBeamData* tpcast = block["E:TPCAST"];
  const RawBeamData* baft1  = block["E:BAFT1"];
  const RawBeamData* baft2  = block["E:BAFT2"];
  const RawBeamData* tortgt = block["E:TORTGT"];
  const RawBeamData *Etrtgtd = block["E:TRTGTD"];

  if(Etrtgtd && Etrtgtd->GetDataLength()) {
    double P2tgt  = Etrtgtd->GetData()[0];
    if(P2tgt < 0.1) {
      //If Beam intensity low, get out.
      return;
    }
  }
  
  if (!(tpcast && baft1 && baft2 && tortgt)) return;
  if (!(tpcast->GetDataLength() &&
	baft1->GetDataLength() &&
	baft2->GetDataLength() &&
	tortgt->GetDataLength()))
    return;
  
  StripHist* sh1 = fStripHist["Baffle Scrape 1"];
  StripHist* sh2 = fStripHist["Baffle Scrape 2"];
  
  double t = tpcast->GetData()[0];
  double pot = tortgt->GetData()[0];
  double y1 = Munits::ToCelcius(Munits::FromFahrenheit(baft1->GetData()[0] - t))/pot;
  double y2 = Munits::ToCelcius(Munits::FromFahrenheit(baft2->GetData()[0] - t))/pot;
  double dae = tortgt->GetSeconds() + 1.0e-6*tortgt->GetMsecs();
  
  sh1->Fill(dae,y1);
  sh2->Fill(dae,y2);

  
  float feet[40] = {-14.76, -2.02, 
		    27.55,   38.55,   49.57,   61.80,  72.80, 
		    83.80,   96.61,   108.84,  133.79, 
		    148.05,  204.60,  261.15,  317.78, 330.13, 
		    351.21,  374.22,  386.34,  407.42, 430.31, 
		    437.55,  458.63,  487.31,  502.01, 519.6, 
		    737.51,  752.20,  814.51,  843.55, 923.94, 
		    938.60,  945.05,  976.62,  997.70, 1020.59, 
		    1027.83, 1048.91, 1070.12, 1081.53};

  std::string lossdevice[40] = {"E:LMV100","E:LMQ101",
				"E:LM1011","E:LM1012","E:LMQ102","E:LM1013","E:LM1014",
				"E:LM1015","E:LMQ103","E:LM1016","E:LMQ104",
				"E:LMQ105","E:LMQ106","E:LMQ107","E:LMQ108","E:LM1081",
				"E:LM1082","E:LMQ109","E:LM1083","E:LM1084","E:LMQ110",
				"E:LM1085","E:LM1086","E:LMQ111","E:LMQ112","E:LMBBBP",
				"E:LMQ113","E:LMQ114","E:LMQ115","E:LMQ116","E:LMQ117",
				"E:LMH117","E:LMQ118","E:LM1181","E:LM1182","E:LMQ119",
				"E:LM1183","E:LM1184","E:LMQ120","E:LMQ121"};
  
  TH1F *sh3 = fTH1F["Beam Loss"];
  sh3->Reset();
  for(int i=0; i<40; i++) {
    const RawBeamData* tmp = block[lossdevice[i]];
    double blah = tmp->GetData()[0];
    sh3->Fill(feet[i], blah );
  }
  
  float feet2[23] = {16.79, 17.70, 60.27, 107.31, 132.58, 145.21, 203.39,
		     258.31, 314.94, 385.04, 436.18,  498.16, 512.88, 736.29,
		     762.90, 811.78, 842.33, 921.10, 943.83, 1095.81, 1096.89,
		     1134.89, 1135.80};
  
  TGraph *sh10 = fTGraph["Beam Position H"];
  TGraph *sh11 = fTGraph["Beam Position V"];
  //MSG("BD",Msg::kDebug) << "A) sh10 has " << sh10->GetN() << "entries" << endl;
  //MSG("BD",Msg::kDebug) << "A) sh11 has " << sh11->GetN() << "entries" << endl;

  for(int i = 0; i < 13; i++) { //sh10->GetN(); i++) {
    sh10->RemovePoint(0);
  }
  for(int j = 0; j < 10; j++) {//sh11->GetN(); j++) {
    sh11->RemovePoint(0);
  }
  
  //MSG("BD",Msg::kDebug) << "B) sh10 has " << sh10->GetN() << "entries" << endl;
  //MSG("BD",Msg::kDebug) << "B) sh11 has " << sh11->GetN() << "entries" << endl;


  ///  List of BPM positions
    std::string bpm_positions[23] = {
      "E:VP101",    "E:HP101",    "E:HP102",    "E:VP103",    "E:HP104",    
      "E:HP105",    "E:VP106",    "E:HP107",    "E:VP108",    "E:HP109",    
      "E:VP110",    "E:VP111",    "E:HP112",    "E:VP113",    "E:HP114", 
      "E:HP115",    "E:VP116",    "E:HP117",    "E:VP118",    "E:HP119",  
      "E:HP121",    "E:HPTGT",    "E:VPTGT"};
    
    for(int i=0; i<23; i++) {
      const RawBeamData* datum = block[bpm_positions[i]];
      if(datum && datum->GetDataLength()) {
	double datum2 = datum->GetData()[0];
	if(datum2 < 990) {
	  /*
	  Explanation for the above cut
	  
	  FYI, according to Peter Lucas:
	  
	  "We now understand the rationale for zero intensity leading to large
	  values of position. This is intentional and is connected with the BPM
	  averaging mask. when there is no beam in the masked-in batches, 1
	  meter is returned as the position. - P. Lucas"
	  */
	  std::string tempor = bpm_positions[i];
	  char ch = tempor[2];
	  //MSG("BD",Msg::kDebug) << "ch=[" << ch << "]" << endl; 
	  if(ch == 'V') {
	    sh11->SetPoint(sh11->GetN(), feet2[i], datum2);
	  } else {
	    sh10->SetPoint(sh10->GetN(), feet2[i], datum2);
	  }
	}
      }
    }

    //Since using TGraph's, have to do this every update.  blah.
    TH1F* frame = sh10->GetHistogram();
    TAxis* axis_x = frame->GetXaxis();
    axis_x->SetTitle("Distance from Q806 (ft)");
    axis_x->CenterTitle();
    TAxis *axis_y = frame->GetYaxis();
    axis_y->SetTitle("Beam Position (mm)");
    axis_y->CenterTitle();
    
    TH1F* frame2 = sh11->GetHistogram();
    TAxis* axis_x2 = frame2->GetXaxis();
    axis_x2->SetTitle("Distance from Q806 (ft)");
    axis_x2->CenterTitle();
    TAxis *axis_y2 = frame2->GetYaxis();
    axis_y2->SetTitle("Beam Position (mm)");
    axis_y2->CenterTitle();

    double min1 = 0.0, max1 = 0.0;
    RangeFinder(sh10, sh11, min1, max1);
    sh10->GetYaxis()->SetRangeUser(min1, max1);
    sh11->GetYaxis()->SetRangeUser(min1, max1);

    const RawBeamData* datum3 = block["I:BEAM"];
    const RawBeamData* datum4 = block["E:TR101D"];
    const RawBeamData* datum5 = block["E:TRTGTD"];

    double datum6 = -1.0, datum7 = -1.0, datum8 = -1.0;

    if(datum3 && datum3->GetDataLength()) {
      datum6 = datum3->GetData()[0];
    }
    if(datum4 && datum4->GetDataLength()) {
      datum7 = datum4->GetData()[0];
    }
    if(datum5 && datum5->GetDataLength()) {
      datum8 = datum5->GetData()[0];
    }
    
    HistMan hm = this->GetHistMan();
    TCanvas *can = (TCanvas*) hm.GetObject("Summary Canvases/Proton Beam Transport Line Summary Canvas");
    TVirtualPad *pad = can->cd(3);
    pad->Clear();
    TLatex l;
    l.SetTextSize(0.1);
    l.DrawLatex(0.1,0.7,Form("Beam From Main Injector = %f x10^{12} (I:BEAM)",datum6));
    l.DrawLatex(0.1,0.4,Form("Beam into NuMI Line = %f x10^{12} (E:TR101D)",datum7));
    l.DrawLatex(0.1,0.1,Form("Beam at NuMI Target = %f x10^{12} (E:TRTGTD)",datum8));
}

void LossModule::RangeFinder(TGraph *sh1, TGraph *sh2, double &min, double &max)
{
  double max1 = GetMax(sh1), max2 = GetMax(sh2);
  double min1 = GetMin(sh1), min2 = GetMin(sh2);

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

  if(min != 0.0) {
    if(min < 0) {
      min *= 1.5;
    } else {
      min *= 0.5;
    }
  } else {
    min = -0.5;
  }

  if(max != 0.0) {
    if(max < 0) {
      max *= 0.5;
    } else {
      max *= 1.5;
    }
  } else {
    max = 0.5;
  }

}

double LossModule::GetMax(TGraph *gr)
{
    if(gr->GetN() == 0) return 0.0;
    double max = 0.0, x = 0.0, y = 0.0;
    for(int i=0; i<gr->GetN(); i++) {
        gr->GetPoint(i, x, y);
        if(y > max) max = y;
    }
    return max;
}
 
double LossModule::GetMin(TGraph *gr)
{
    if(gr->GetN() == 0) return 0.0;
    double max = GetMax(gr);
    double min = max, x = 0.0, y = 0.0;
    for(int i=0; i<gr->GetN(); i++) {
        gr->GetPoint(i, x, y);
        if(y < min) min = y;
    }
    return min;
}