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

#include "HadMuMonModule.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/RawBeamMonHeaderBlock.h>
#include <RawData/RawBeamData.h>

#include <BeamDataUtil/BDSwicCalibrator.h>
#include <BeamDataUtil/BDHadMuMon.h>
#include <BeamDataUtil/BDDevices.h>
#include <BeamDataUtil/BDSwicDevice.h>

#include <DataUtil/GetRawBlock.h>
#include <Conventions/Munits.h>

#include <TPad.h>
#include <TH2F.h>
#include <TGraph.h>
#include <TMultiGraph.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TPaletteAxis.h>

//ClassImp(HadMuMonModule)

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

using namespace std;
using namespace DataUtil;

CVSID("$Id: HadMuMonModule.cxx,v 1.9 2005/06/09 15:26:43 thosieck Exp $");
JOBMODULE(HadMuMonModule,"HadMuMon","HadMuMon quantities related for Monitoring");

HadMuMonModule::HadMuMonModule()
{

  StripHist *sh1  = new StripHist("Beam Intensity Alcove 1 (mV)","Muon Monitor Beam Intensity Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh2  = new StripHist("Beam Intensity Alcove 2 (mV)","Muon Monitor Beam Intensity Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh3  = new StripHist("Beam Intensity Alcove 3 (mV)","Muon Monitor Beam Intensity Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh4  = new StripHist("Beam X Centroid Alcove 1","Beam X Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh5  = new StripHist("Beam X Centroid Alcove 2","Beam X Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh6  = new StripHist("Beam X Centroid Alcove 3","Beam X Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh7  = new StripHist("Beam Y Centroid Alcove 1","Beam Y Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh8  = new StripHist("Beam Y Centroid Alcove 2","Beam Y Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh9  = new StripHist("Beam Y Centroid Alcove 3","Beam Y Centroid Red=Alcove 1 Blue=Alcove 2 Black=Alcove 3");
  StripHist *sh10 = new StripHist("Beam Intensity HadMon (mV)","Hadron Monitor Beam Intensity");
  StripHist *sh11 = new StripHist("Beam X HadMon","Hadron Monitor Beam Position Red=X Blue=Y");
  StripHist *sh12 = new StripHist("Beam Y HadMon","Hadron Monitor Beam Position Red=X Blue=Y");
  StripHist *sh13 = new StripHist("Beam X RMS HadMon","Hadron Monitor Beam RMS Red=X Blue=Y");
  StripHist *sh14 = new StripHist("Beam Y RMS HadMon","Hadron Monitor Beam RMS Red=X Blue=Y");
  
  TH2F *MuonMonitor1  = new TH2F("MuonMonitor1","Alcove 1 Muon Monitor Intensity",9,-45,45,9,-45,45);
  TH2F *MuonMonitor2  = new TH2F("MuonMonitor2","Alcove 2 Muon Monitor Intensity",9,-45,45,9,-45,45);
  TH2F *MuonMonitor3  = new TH2F("MuonMonitor3","Alcove 3 Muon Monitor Intensity",9,-45,45,9,-45,45);
  TH2F *HadronMonitor = new TH2F("HadronMonitor","Hadron Monitor Intensity",7,-15.75,15.75,7,-15.75,15.75);

  TH1F *HadronMonitorX = new TH1F("HadronMonitorX","Hadron Monitor Intensity in X Position",7,-15.75,15.75);
  TH1F *HadronMonitorY = new TH1F("HadronMonitorY","Hadron Monitor Intensity in Y Position",7,-15.75,15.75);
  
  sh1  -> SetLineColor(2);
  sh2  -> SetLineColor(4);
  sh4  -> SetLineColor(2);
  sh5  -> SetLineColor(4);
  sh7  -> SetLineColor(2);
  sh8  -> SetLineColor(4);
  sh11 -> SetLineColor(2);
  sh12 -> SetLineColor(4);
  sh13 -> SetLineColor(2);
  sh14 -> SetLineColor(4);

  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);
  sh7  -> SetStripRange(1*Munits::day);
  sh8  -> SetStripRange(1*Munits::day);
  sh9  -> SetStripRange(1*Munits::day);
  sh10 -> SetStripRange(1*Munits::day);
  sh11 -> SetStripRange(1*Munits::day);
  sh12 -> SetStripRange(1*Munits::day);
  sh13 -> SetStripRange(1*Munits::day);
  sh14 -> SetStripRange(1*Munits::day);
  
  HadronMonitorX->GetXaxis()->SetTitle("Horizontal Position (inches)");
  HadronMonitorY->GetXaxis()->SetTitle("Vertical Position (inches)");

  HadronMonitorX->GetYaxis()->SetTitle("Beam Intensity (mV)");
  HadronMonitorY->GetYaxis()->SetTitle("Beam Intensity (mV)");

  HadronMonitor->GetXaxis()->SetTitle("Horizontal Position (inches)");
  HadronMonitor->GetYaxis()->SetTitle("Vertical Position (inches)");

  MuonMonitor1->GetXaxis()->SetTitle("Horizontal Position (inches)");
  MuonMonitor1->GetYaxis()->SetTitle("Vertical Position (inches)");
  
  MuonMonitor2->GetXaxis()->SetTitle("Horizontal Position (inches)");
  MuonMonitor2->GetYaxis()->SetTitle("Vertical Position (inches)");

  MuonMonitor3->GetXaxis()->SetTitle("Horizontal Position (inches)");
  MuonMonitor3->GetYaxis()->SetTitle("Vertical Position (inches)");

  HadronMonitor->GetXaxis()->CenterTitle();
  HadronMonitor->GetYaxis()->CenterTitle();
  
  MuonMonitor1->GetXaxis()->CenterTitle();
  MuonMonitor1->GetYaxis()->CenterTitle();

  MuonMonitor2->GetXaxis()->CenterTitle();
  MuonMonitor2->GetYaxis()->CenterTitle();

  MuonMonitor3->GetXaxis()->CenterTitle();
  MuonMonitor3->GetYaxis()->CenterTitle();

  HadronMonitorX->GetXaxis()->CenterTitle();
  HadronMonitorY->GetXaxis()->CenterTitle();

  HadronMonitorX->GetYaxis()->CenterTitle();
  HadronMonitorY->GetYaxis()->CenterTitle();

  MuonMonitor1->SetStats(0);
  MuonMonitor2->SetStats(0);
  MuonMonitor3->SetStats(0);
  HadronMonitor->SetStats(0);

  sh1->SetGraphYTitle("Beam Intensity (mV)");
  sh2->SetGraphYTitle("Beam Intensity (mV)");
  sh3->SetGraphYTitle("Beam Intensity (mV)");

  sh4->SetGraphYTitle("Beam X Centroid (mm)");
  sh5->SetGraphYTitle("Beam X Centroid (mm)");
  sh6->SetGraphYTitle("Beam X Centroid (mm)");
  
  sh7->SetGraphYTitle("Beam Y Centroid (mm)");
  sh8->SetGraphYTitle("Beam Y Centroid (mm)");
  sh9->SetGraphYTitle("Beam Y Centroid (mm)");

  sh10->SetGraphYTitle("Beam Intensity (mV)");
  sh11->SetGraphYTitle("Beam Position (mm)");
  sh12->SetGraphYTitle("Beam Position (mm)");
  sh13->SetGraphYTitle("Beam RMS (mm)");
  sh14->SetGraphYTitle("Beam RMS (mm)");

  fStripHist["Beam Intensity Alcove 1"]  = sh1;
  fStripHist["Beam Intensity Alcove 2"]  = sh2;
  fStripHist["Beam Intensity Alcove 3"]  = sh3;
  fStripHist["Beam X Centroid Alcove 1"] = sh4;
  fStripHist["Beam X Centroid Alcove 2"] = sh5;
  fStripHist["Beam X Centroid Alcove 3"] = sh6;
  fStripHist["Beam Y Centroid Alcove 1"] = sh7;
  fStripHist["Beam Y Centroid Alcove 2"] = sh8;
  fStripHist["Beam Y Centroid Alcove 3"] = sh9;  
  fStripHist["Beam Intensity HadMon"]    = sh10;
  fStripHist["Beam X HadMon"]            = sh11;
  fStripHist["Beam Y HadMon"]            = sh12;
  fStripHist["Beam X RMS HadMon"]        = sh13;
  fStripHist["Beam Y RMS HadMon"]        = sh14;

  fTH2F["MuonMonitor1"]  = MuonMonitor1;
  fTH2F["MuonMonitor2"]  = MuonMonitor2;
  fTH2F["MuonMonitor3"]  = MuonMonitor3;
  fTH2F["HadronMonitor"] = HadronMonitor;

  fTH1F["HadronMonitorX"] = HadronMonitorX;
  fTH1F["HadronMonitorY"] = HadronMonitorY;
}

HadMuMonModule::~HadMuMonModule()
{

}

typedef std::map<std::string, TH2F*> TH2FMap;
typedef std::map<std::string, TH1F*> TH1FMap;

void HadMuMonModule::BeginJob()
{
  MSG("BD",Msg::kDebug) << "Adding Histograms!" << endl;

  HistMan hm = this->GetHistMan();
  
  StripHist *sh1 = fStripHist["Beam Intensity Alcove 1"];
  StripHist *sh2 = fStripHist["Beam Intensity Alcove 2"];
  StripHist *sh3 = fStripHist["Beam Intensity Alcove 3"];

  TCanvas *canvas = new TCanvas("Beam Intensity from Muon Monitors","Muon Monitor Intensity",500,400);

  sh1->DrawStrip("AL");
  sh2->DrawStrip("LP");
  sh3->DrawStrip("L");
  
  hm.Adopt("HadMuMon",canvas);

  StripHist *sh4 = fStripHist["Beam X Centroid Alcove 1"];
  StripHist *sh5 = fStripHist["Beam X Centroid Alcove 2"];
  StripHist *sh6 = fStripHist["Beam X Centroid Alcove 3"];
  
  TCanvas *canvas2 = new TCanvas("Beam X Centroid from Muon Monitors","X Centroid",500,400);

  sh4->DrawStrip("AL");
  sh5->DrawStrip("LP");
  sh6->DrawStrip("L");
  
  hm.Adopt("HadMuMon",canvas2);

  StripHist *sh7 = fStripHist["Beam Y Centroid Alcove 1"];
  StripHist *sh8 = fStripHist["Beam Y Centroid Alcove 2"];
  StripHist *sh9 = fStripHist["Beam Y Centroid Alcove 3"];
  
  TCanvas *canvas3 = new TCanvas("Beam Y Centroid from Muon Monitors","Y Centroid",500,400);

  sh7->DrawStrip("AL");
  sh8->DrawStrip("LP");
  sh9->DrawStrip("L");
  
  hm.Adopt("HadMuMon",canvas3);

  StripHist *sh10 = fStripHist["Beam Intensity HadMon"];

  TCanvas *canvas4 = new TCanvas("Beam Intensity from Hadron Monitor","Hadron Monitor Intensity",500,400);

  sh10->DrawStrip("AL");

  hm.Adopt("HadMuMon",canvas4);

  StripHist *sh11 = fStripHist["Beam X HadMon"];
  StripHist *sh12 = fStripHist["Beam Y HadMon"];
  
  TCanvas *canvas5 = new TCanvas("Hadron Monitor Beam Location","Beam Location",500,400);
  
  sh11->DrawStrip("AL");
  sh12->DrawStrip("LP");
  
  hm.Adopt("HadMuMon",canvas5);

  StripHist *sh13 = fStripHist["Beam X RMS HadMon"];
  StripHist *sh14 = fStripHist["Beam Y RMS HadMon"];
  
  TCanvas *canvas6 = new TCanvas("Hadron Monitor Beam Size","Beam Size",500,400);
  
  sh13->DrawStrip("AL");
  sh14->DrawStrip("LP");

  hm.Adopt("HadMuMon",canvas6);
  
  TH2FMap::iterator uchicago, done2 = fTH2F.end();
  for( uchicago = fTH2F.begin(); uchicago != done2; ++ uchicago) {
    TH2F *hist = uchicago->second;
    TCanvas* canvas = new TCanvas(hist->GetTitle(),hist->GetTitle(),500,400);
    hist->Draw("colz");
    hm.Adopt("HadMuMon",canvas);
  }
  
  TH1FMap::iterator caltech, done3 = fTH1F.end();
  for( caltech = fTH1F.begin(); caltech != done3; ++caltech) {
    TH1F *hist = caltech->second;
    TCanvas* canvas8 = new TCanvas(hist->GetTitle(),hist->GetTitle(),500,400);
    hist->Draw();
    hm.Adopt("HadMuMon",canvas8);
  }  

  TH2F *th1 = fTH2F["MuonMonitor1"];
  TH2F *th2 = fTH2F["MuonMonitor2"];
  TH2F *th3 = fTH2F["MuonMonitor3"];
  TH2F *th4 = fTH2F["HadronMonitor"];
  TH1F *th5 = fTH1F["HadronMonitorX"];
  TH1F *th6 = fTH1F["HadronMonitorY"];

  TCanvas *MuMonSumCan = new TCanvas("Muon Monitor Summary Canvas","Muon Monitor Summary Canvas",800,800);
  MuMonSumCan -> Divide(2,3);
  MuMonSumCan -> cd(1);
  th1 -> Draw("colz");
  MuMonSumCan -> cd(2);
  sh1 -> DrawStrip("AL");
  sh2 -> DrawStrip("LP");
  sh3 -> DrawStrip("L");
  MuMonSumCan -> cd(3);
  th2 -> Draw("colz");
  MuMonSumCan -> cd(4);
  sh4 -> DrawStrip("AL");
  sh5 -> DrawStrip("LP");
  sh6 -> DrawStrip("L");
  MuMonSumCan -> cd(5);
  th3 -> Draw("colz");
  MuMonSumCan -> cd(6);
  sh7 -> DrawStrip("AL");
  sh8 -> DrawStrip("LP");
  sh9 -> DrawStrip("L");

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

  TCanvas *HadMonSumCan = new TCanvas("Hadron Monitor Summary Canvas","Hadron Monitor Summary Canvas",800,800);
  HadMonSumCan -> Divide(2,3);
  HadMonSumCan -> cd(1);
  th4->Draw("colz");
  HadMonSumCan -> cd(2);
  sh10->DrawStrip("AL");
  HadMonSumCan -> cd(3);
  th5->Draw();
  HadMonSumCan -> cd(4);
  sh11->DrawStrip("AL");
  sh12->DrawStrip("L");
  HadMonSumCan -> cd(5);
  th6->Draw();
  HadMonSumCan -> cd(6);
  sh13->DrawStrip("AL");
  sh14->DrawStrip("L");

  hm.Adopt("Summary Canvases",HadMonSumCan);
   
}

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

  const RawBeamData *mma1ds = block["E:MMA1DS"]; 
  const RawBeamData *mma2ds = block["E:MMA2DS"]; 
  const RawBeamData *mma3ds = block["E:MMA3DS"];
  const RawBeamData *hadmds = block["E:HADMDS"];
  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(!(mma1ds && mma2ds && mma3ds)) {
    MSG("BD",Msg::kDebug) << "Do not have Muon Monitor Data!" << endl;
    if(!mma1ds) MSG("BD",Msg::kDebug) << "  No MMA1DS" << endl;
    if(!mma2ds) MSG("BD",Msg::kDebug) << "  No MMA2DS" << endl;
    if(!mma3ds) MSG("BD",Msg::kDebug) << "  No MMA3DS" << endl;
    return;
  }

  if(!(mma1ds->GetDataLength() &&
       mma2ds->GetDataLength() &&
       mma3ds->GetDataLength())) {
    MSG("BD",Msg::kDebug) << "  Muon Monitor device but no data!" << endl;
    return;
  }

  if(!hadmds || !hadmds->GetDataLength()) {
    MSG("BD",Msg::kDebug) << "  No Hadron Monitor Data!" << endl;
    return;
  }

  MSG("BD",Msg::kDebug) << mma1ds->GetName() << ", len = " << mma1ds->GetDataLength() << endl;
  BDHadMuMon hmm1(*mma1ds);
  MSG("BD",Msg::kDebug) << mma1ds->GetName() << ", len = " << mma1ds->GetDataLength() << endl;
  MSG("BD",Msg::kDebug) << mma1ds->GetName() << ", tot = " << hmm1.GetTotalVoltage() << endl;
  BDHadMuMon hmm2(*mma2ds);
  MSG("BD",Msg::kDebug) << mma2ds->GetName() << ", tot = " << hmm2.GetTotalVoltage() << endl;
  BDHadMuMon hmm3(*mma3ds);
  MSG("BD",Msg::kDebug) << mma3ds->GetName() << ", tot = " << hmm3.GetTotalVoltage() << endl;
  BDHadMuMon hmm4(*hadmds);
  MSG("BD",Msg::kDebug) << hadmds->GetName() << ", tot = " << hmm4.GetTotalVoltage() << endl;

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

  TH2F *h1 = fTH2F["MuonMonitor1"];
  TH2F *h2 = fTH2F["MuonMonitor2"];
  TH2F *h3 = fTH2F["MuonMonitor3"];
  TH2F *h4 = fTH2F["HadronMonitor"];

  h1->Reset(); h2->Reset(); h3->Reset(); h4->Reset();

  TH1F *h5 = fTH1F["HadronMonitorX"];
  TH1F *h6 = fTH1F["HadronMonitorY"];

  h5->Reset(); h6->Reset();

  StripHist* sh1  = fStripHist["Beam Intensity Alcove 1"];
  StripHist* sh2  = fStripHist["Beam Intensity Alcove 2"];
  StripHist* sh3  = fStripHist["Beam Intensity Alcove 3"];
  StripHist* sh4  = fStripHist["Beam X Centroid Alcove 1"];
  StripHist* sh5  = fStripHist["Beam X Centroid Alcove 2"];
  StripHist* sh6  = fStripHist["Beam X Centroid Alcove 3"];
  StripHist* sh7  = fStripHist["Beam Y Centroid Alcove 1"];
  StripHist* sh8  = fStripHist["Beam Y Centroid Alcove 2"];
  StripHist* sh9  = fStripHist["Beam Y Centroid Alcove 3"];
  StripHist* sh10 = fStripHist["Beam Intensity HadMon"];
  StripHist* sh11 = fStripHist["Beam X HadMon"];
  StripHist* sh12 = fStripHist["Beam Y HadMon"];
  StripHist* sh13 = fStripHist["Beam X RMS HadMon"];
  StripHist* sh14 = fStripHist["Beam Y RMS HadMon"];
  
  for(int row = 1; row <= hmm1.GetNrowcol(); row++) {
    for(int col = 1; col <= hmm1.GetNrowcol(); col++) {
      // muon monitors all have identical channel mappings
      int index = hmm1.Index(hmm1.Channel(row,col));
      h1->Fill(hmm1.PixelPosition(col) / Munits::inch, 
	       hmm1.PixelPosition(row) / Munits::inch, 
	       hmm1.GetVoltage(index)/Munits::millivolt);
      
      h2->Fill(hmm2.PixelPosition(col) / Munits::inch, 
	       hmm2.PixelPosition(row) / Munits::inch, 
	       hmm2.GetVoltage(index)/Munits::millivolt);

      h3->Fill(hmm3.PixelPosition(col) / Munits::inch, 
	       hmm3.PixelPosition(row) / Munits::inch, 
	       hmm3.GetVoltage(index)/Munits::millivolt);
    }
  }

  for(int row = 1; row <= hmm4.GetNrowcol(); row++) {
    for(int col = 1; col <= hmm4.GetNrowcol(); col++) {
      double v = hmm4.GetVoltage(hmm4.Index(hmm4.Channel(row, col)));
      v /= Munits::millivolt;
      h4->Fill(hmm4.PixelPosition(col) / Munits::inch,
	       hmm4.PixelPosition(row) / Munits::inch, v);
      h5->Fill(hmm4.PixelPosition(col) / Munits::inch, v);
      h6->Fill(hmm4.PixelPosition(row) / Munits::inch, v);
    }
  }

  h5->GetYaxis()->SetTitleOffset(1.2);
  h6->GetYaxis()->SetTitleOffset(1.2);

  double q1   = hmm1.GetTotalVoltage() / Munits::millivolt;
  double q2   = hmm2.GetTotalVoltage() / Munits::millivolt;
  double q3   = hmm3.GetTotalVoltage() / Munits::millivolt;
  double q4   = hmm4.GetTotalVoltage() / Munits::millivolt;
  double dae1 = mma1ds->GetSeconds() + 1.0e-6*mma1ds->GetMsecs();
  double dae2 = mma2ds->GetSeconds() + 1.0e-6*mma2ds->GetMsecs();
  double dae3 = mma3ds->GetSeconds() + 1.0e-6*mma3ds->GetMsecs();
  double dae4 = hadmds->GetSeconds() + 1.0e-6*hadmds->GetMsecs();

  sh1  -> Fill(dae1, q1);
  sh2  -> Fill(dae2, q2);
  sh3  -> Fill(dae3, q3);
  sh10 -> Fill(dae4, q4);

  double min1 = 0.0, max1 = 1.0;

  RangeFinder(sh1, sh2, sh3, min1, max1);

  sh1->GetStrip().GetYaxis()->SetRangeUser(min1, max1);
  sh2->GetStrip().GetYaxis()->SetRangeUser(min1, max1);
  sh3->GetStrip().GetYaxis()->SetRangeUser(min1, max1);
  
  double xmean1 = 0.0, ymean1 = 0.0, xrms1 = 0.0, yrms1 = 0.0;
  double xmean2 = 0.0, ymean2 = 0.0, xrms2 = 0.0, yrms2 = 0.0;
  double xmean3 = 0.0, ymean3 = 0.0, xrms3 = 0.0, yrms3 = 0.0;
  double xmean4 = 0.0, ymean4 = 0.0, xrms4 = 0.0, yrms4 = 0.0;
  
  hmm1.GetStats(xmean1, ymean1, xrms1, yrms1);
  hmm2.GetStats(xmean2, ymean2, xrms2, yrms2);
  hmm3.GetStats(xmean3, ymean3, xrms3, yrms3);
  hmm4.GetStats(xmean4, ymean4, xrms4, yrms4);
  
  sh4 -> Fill(dae1, xmean1 / Munits::mm);
  sh5 -> Fill(dae2, xmean2 / Munits::mm);
  sh6 -> Fill(dae3, xmean3 / Munits::mm);

  double min2 = 0.0, max2 = 1.0;

  RangeFinder(sh4, sh5, sh6, min2, max2);

  sh4->GetStrip().GetYaxis()->SetRangeUser(min2, max2);
  sh5->GetStrip().GetYaxis()->SetRangeUser(min2, max2);
  sh6->GetStrip().GetYaxis()->SetRangeUser(min2, max2);

  sh7  -> Fill(dae1, ymean1 / Munits::mm);
  sh8  -> Fill(dae2, ymean2 / Munits::mm);
  sh9  -> Fill(dae3, ymean3 / Munits::mm);

  double min3 = 0.0, max3 = 1.0;

  RangeFinder(sh7, sh8, sh9, min3, max3);

  sh7->GetStrip().GetYaxis()->SetRangeUser(min3, max3);
  sh8->GetStrip().GetYaxis()->SetRangeUser(min3, max3);
  sh9->GetStrip().GetYaxis()->SetRangeUser(min3, max3);

  sh11 -> Fill(dae4, xmean4 / Munits::mm);
  sh12 -> Fill(dae4, ymean4 / Munits::mm);

  double min4 = 0.0, max4 = 1.0;

  RangeFinder(sh11, sh12, min4, max4);

  sh11->GetStrip().GetYaxis()->SetRangeUser(min4, max4);
  sh12->GetStrip().GetYaxis()->SetRangeUser(min4, max4);

  sh13 -> Fill(dae4, xrms4  / Munits::mm);
  sh14 -> Fill(dae4, yrms4  / Munits::mm);

  double min5 = 0.0, max5 = 1.0;

  RangeFinder(sh13, sh14, min5, max5);

  sh13->GetStrip().GetYaxis()->SetRangeUser(min5, max5);
  sh14->GetStrip().GetYaxis()->SetRangeUser(min5, max5);

}

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

  if(min1 < min2) {
    if(min1 < min3) {
      min = min1;
    } else {
      min = min3;
    }
  } else {
    if(min2 < min3) {
      min = min2;
    } else {
      min = min3;
    }
  }

  if(max1 > max2) {
    if(max1 > max3) {
      max = max1;
    } else {
      max = max3;
    }
  } else {
    if(max2 > max3) {
      max = max2;
    } else {
      max = max3;
    }
  }

  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;
  }

}

void HadMuMonModule::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(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;
  }

}