NCDataQualityModule.cxx

Go to the documentation of this file.

#include "NCDataQualityModule.h"

#include "TArc.h"
#include "TArrow.h"
#include "TCanvas.h"
#include "TF1.h"
#include "TFile.h"
#include "TGaxis.h"
#include "TGraph.h"
#include "TGraphAsymmErrors.h"
#include "TGraphErrors.h"
#include "TLegend.h"
#include "TLatex.h"
#include "TList.h"
#include "TObjString.h"
#include "TProfile.h"
#include "TSystem.h"
#include "TSystemDirectory.h"

#include "MessageService/MsgService.h"

#include "JobControl/JobCModuleRegistry.h" // For JOBMODULE macro

#include "NCUtils/NCRunUtil.h"
#include "NCUtils/NCUtility.h"

#include <cassert>
#include <iostream>
using namespace NCType;
using std::vector;
using NC::Utility::SQR;

CVSID("$Id: NCDataQualityModule.cxx,v 1.20 2009/07/28 12:10:04 bckhouse Exp $");


// Declare this module to JobControl. Arguments are:
//  (1) The class name
//  (2) The human-readable name
//  (3) A short, human-readable description of what the module does
JOBMODULE(NCDataQualityModule,
          "NCDataQualityModule",
          "Get DQ plots out of files and draw them");


NCDataQualityModule::NCDataQualityModule():
  fNearPOTTotal(-1), fFarPOTTotal(-1),
  fNMCFilesNear(0), fNMCFilesFar(0)

{
  for(int n = 0; n < kDQNumCombinations; ++n) fDQPlots[n] = 0;

  for(int k = 0; k < kDQ2DNumDists; ++k){
    fDQNearData2DTotal.push_back(new NCDataQualityPlot2D(k, "NearDataTotal"));
    fDQNearData2DNC.push_back(new NCDataQualityPlot2D(k, "NearDataNC"));
    fDQNearData2DCC.push_back(new NCDataQualityPlot2D(k, "NearDataCC"));
  }
}

NCDataQualityModule::~NCDataQualityModule()
{
  for(int n = 0; n < kDQNumCombinations; ++n)
    if(fDQPlots[n]) delete fDQPlots[n];
}

//......................................................................
const Registry& NCDataQualityModule::DefaultConfig() const
{
  static Registry r;

  r.UnLockValues();

  r.Set("DataMCPath",          "dataFiles/*.root");
  r.Set("UseMCAsData",         false);
  r.Set("UseCCEvents",         true);
  r.Set("UseNCEvents",         true);
  r.Set("UseMockData",         false);
  r.Set("MockDataSet",         "F21910001");

  TString runsToUse;
  for(int n = 0; n <= NC::RunUtil::kMaxRun; ++n)
    runsToUse += TString::Format("%d ", n);
  r.Set("RunsToUse",              runsToUse);

  r.Set("BeamType",            "L010z185i");
  r.Set("FileLimitNearData",   1000000000);
  r.Set("FileLimitFarData",    10000);
  r.Set("FileLimitNearMC",     1000000000);
  r.Set("FileLimitFarMC",      10000);
  r.Set("MockDataSubRun",      0);
  r.Set("ScaleNearExposures",  false);
  r.Set("FileName",            "dqplots.root");

  r.LockValues();
  return r;
}


//......................................................................
void NCDataQualityModule::Config(const Registry& r)
{
  int         tmpb;
  int         tmpi;
  const char* tmps;

  if (r.Get("DataMCPath",            tmps)){
    //for some reason the TSystemDirectory constructor cant resolve
    //environmental variables, so just take care of it here
    int pos = 0;
    TString temp = tmps;
    if(temp.Contains("$")){
      pos = temp.Index("/");
      temp.Resize(pos);
      temp.Remove(TString::kLeading, '$');
      fDataMCPath                  = tmps;
      fDataMCPath.Replace(0, pos, gSystem->Getenv(temp));
    }
    else fDataMCPath = temp;
  }

  if(r.Get("UseMCAsData",           tmpb)) fUseMCAsData         = tmpb;
  if(r.Get("UseCCEvents",           tmpb)) fUseCC               = tmpb;
  if(r.Get("UseNCEvents",           tmpb)) fUseNC               = tmpb;
  if(r.Get("UseMockData",           tmpb)) fUseMockData         = tmpb;
  if(r.Get("MockDataSet",           tmps)) fMockDataSet         = tmps;

  if(r.Get("RunsToUse",             tmps)){
    fRunsToUse.clear();
    vector<int> runs = NC::Utility::ParseNumberList(tmps);
    for(unsigned int n = 0; n < runs.size(); ++n)
      fRunsToUse.push_back(NC::RunUtil::ERunType(runs[n]));
  }

  if(r.Get("FileLimitNearData",     tmpi)) fFileLimitNearData   = tmpi;
  if(r.Get("FileLimitFarData",      tmpi)) fFileLimitFarData    = tmpi;
  if(r.Get("FileLimitNearMC",       tmpi)) fFileLimitNearMC     = tmpi;
  if(r.Get("FileLimitFarMC",        tmpi)) fFileLimitFarMC      = tmpi;
  if(r.Get("MockDataSubRun",        tmpi)) fMockDataSubRun      = tmpi;
  if(r.Get("ScaleNearExposures",    tmpb)) fScaleNearExposures  = tmpb;
  if(r.Get("FileName",              tmps)) fFileName            = tmps;

  if (r.Get("BeamType",              tmps))
    fBeamIndex = NCType::BeamListFromString(tmps);
}


void NCDataQualityModule::EndJob()
{
  GetConfig().SetName("config_registry");
  GetConfig().Print();

  CombineDataQualityPlots(Detector::kNear, kDQData);
  CombineDataQualityPlots(Detector::kNear, kDQMC);
  CombineDataQualityPlots(Detector::kFar,  kDQData);
  CombineDataQualityPlots(Detector::kFar,  kDQMC);
  CombineDataQualityPlots(Detector::kFar,  kDQMCOsc);

  SetPOTValues();

  TFile f(fFileName, "RECREATE");
  GetConfig().Write();
  DrawDataQualityPlots();
  f.Write();
  f.Close();

  MSG("NCDataQualityModule", Msg::kInfo) << "Output filename "
                                          << fFileName
                                          << endl;
}

void NCDataQualityModule::DrawDataQualityPlots()
{
  // get a pointer to the current directory
  // this is one of the output files
  TDirectory* saveDir = gDirectory;

  //loop over the data quality histograms
  //make folders for data/mc and subfolders for near/far
  TDirectory* dir = saveDir->mkdir("DataQuality", "DataQuality");
  dir->cd();

  TDirectory* near = dir->mkdir("NearDetector", "NearDetector");
  TDirectory* far = dir->mkdir("FarDetector", "FarDetector");

  near->cd();
  DrawNearDataQualityPlots();

  far->cd();
  DrawFarDataQualityPlots();

  saveDir->cd();
}


void NCDataQualityModule::SetPOTValues()
{
  NCPOTCounter potcount;

  potcount.Config(GetConfig(), fBeamIndex, fRunsToUse);

  map<NCBeam::Info, double> dataPOTNear;
  map<NCBeam::Info, double> dataPOTFar;

  const NCType::EFileType fileType = fUseMockData ? NCType::kMockFile : NCType::kBeamFile;

  potcount.SetPOTValues(dataPOTNear, Detector::kNear,
                        NCType::kBeamFile, NCType::kData);
  potcount.SetPOTValues(dataPOTFar, Detector::kFar,
                        fileType, NCType::kData);

  const NCBeam::Info runI (BeamType::kL010z185i, NC::RunUtil::kRunI);
  const NCBeam::Info runII(BeamType::kL010z185i, NC::RunUtil::kRunII);

  fFarPOTTotal = dataPOTFar[runI] + dataPOTFar[runII];
  fNearPOTTotal = dataPOTNear[runI] + dataPOTNear[runII];
}


NCDataQualityPlot* NCDataQualityModule::GetDQPlot(unsigned int dist,
                                                  unsigned int nccc,
                                                  unsigned int datamc,
                                                  Detector::Detector_t det,
                                                  unsigned int intensity)
{
  // No-one should be trying to make ND oscillated MC.
  if(datamc == kDQMCOsc && det == Detector::kNear) return 0;

  assert(nccc < kDQNumInteractions);
  assert(datamc < kDQNumDatasets);
  assert(det == Detector::kNear || det == Detector::kFar);
  assert(intensity < kDQNumIntensities);

  // We are reliant on the actual numerical values of these, so check first.
  assert(Detector::kNear == 1 && Detector::kFar == 2);
  const unsigned int detOffset = det - 1;

  const unsigned int numDetectors = 2;

  const unsigned int totOffset = nccc +
                         kDQNumInteractions*datamc +
                         kDQNumInteractions*kDQNumDatasets*detOffset +
                         kDQNumInteractions*kDQNumDatasets*numDetectors*intensity +
                         kDQNumInteractions*kDQNumDatasets*numDetectors*kDQNumIntensities*dist;

  assert(totOffset < UInt_t(kDQNumCombinations));

  if(!fDQPlots[totOffset])
    fDQPlots[totOffset] = CreateDQPlot(dist, nccc, datamc, det, intensity);

  return fDQPlots[totOffset];
}

NCDataQualityPlot* NCDataQualityModule::CreateDQPlot(unsigned int dist,
                                                     unsigned int nccc,
                                                     unsigned int datamc,
                                                     Detector::Detector_t det,
                                                     unsigned int intense)
{
  TString suffix = Detector::AsString(det);
  suffix += (datamc == kDQData) ? "Data" : "MC";
  if(nccc == kDQNC) suffix += "NC";
  if(nccc == kDQCC) suffix += "CC";
  if(nccc == kDQTotal) suffix += "Total";

  if(intense == kDQLowIntensity) suffix += "LowIntensity";
  if(datamc == kDQMCOsc) suffix += "Osc";

  NCDataQualityPlot* ret = new NCDataQualityPlot(dist, suffix);
  // Otherwise we end up owned by someone who goes out of scope -> bad things
  ret->SetDirectory(0);
  return ret;
}

//-----------------------------------------------------------------------
void NCDataQualityModule::CombineDataQualityPlots(Detector::Detector_t det,
                                                  unsigned int datamc)
{
  // Need a special suffix for oscillated MC
  TString suffix = (datamc == kDQMCOsc) ? "Osc" : "";

  static bool firstTimeNear=true;
  static bool firstTimeFar=true;

  TString fileType = GetFileTypeString(det, datamc);

  MSG("NCDataQualityModule", Msg::kInfo) << "fDataMCPath = " << fDataMCPath << endl;

  TString detName = (det == Detector::kFar) ? "far" : "near";

  TString totalName = "MCTotal"+suffix;
  TString ncName = "MCNC"+suffix;
  TString ccName = "MCCC"+suffix;
  TString name = "";

  //Now loop over the BeamType::AsString values and fill the
  //data and mc chains
  typedef std::vector<BeamType::BeamType_t>::iterator it_t;
  for(it_t beamIter = fBeamIndex.begin(); beamIter != fBeamIndex.end(); ++beamIter){
    BeamType::BeamType_t bType = *beamIter;
    TString beamType = BeamType::AsString(bType);

    //the files for different beam configurations are in uDST/LXXXzYYYi
    //directories.  that is directories named for the beam type
    TSystemDirectory dataDir(beamType, fDataMCPath+beamType+"/");

    MSG("NCDataQualityModule", Msg::kInfo) << "Adding DQ plots from "
                                            << (fDataMCPath+beamType)
                                            << endl;

    TList* files = dataDir.GetListOfFiles();
    assert(files);

    //Loop over files in dataDir
    for(int de = 0; de < files->GetEntries(); ++de) {
      //Do nothing for "." and ".." entries
      TString fileName = files->At(de)->GetName();
      if(fileName == "." || fileName == "..") continue;

      if(fileName.Contains(beamType)
         && fileName.Contains("strip")
         && fileName.Contains(detName)
         && fileName.Contains(fileType)){

        if (datamc==kDQMC) {
          if (det==Detector::kNear && firstTimeNear) ++fNMCFilesNear;
          if (det==Detector::kFar && firstTimeFar) ++fNMCFilesFar;
        }
        TFile tf(fileName);

        const NC::RunUtil::ERunType runType = NC::RunUtil::FindRunType(fileName);
        if(fileType.Contains("data") || fileType.Contains("blind")){
          totalName = "DataTotal";
          ncName = "DataNC";
          ccName = "DataCC";

          bool inRunsToUse = false;
          for(unsigned int n = 0; n < fRunsToUse.size(); ++n)
            if(fRunsToUse[n] == runType) inRunsToUse = true;
          if(!inRunsToUse) continue;

          TGraph* graph = (TGraph*)(tf.Get("DataQuality/vertexXY"));
          if(det == Detector::kFar && graph){
            for(int p = 0; p < graph->GetN(); ++p){
              double x, y;
              graph->GetPoint(p, x, y);
              fFarVtxX.push_back(x);
              fFarVtxY.push_back(y);
            }
          }
          if(det == Detector::kNear && graph){
            for(int p = 0; p < graph->GetN(); ++p){
              double x, y;
              graph->GetPoint(p, x, y);
              fNearVtxX.push_back(x);
              fNearVtxY.push_back(y);
            }
          }
          for(int i = 0; i < kDQ2DNumDists && det == Detector::kNear; ++i){
            NCDataQualityPlot2D* hist2D;
            hist2D = (NCDataQualityPlot2D*)(tf.Get("DataQuality/"+kDQVars2D[i].name+totalName));
            fDQNearData2DTotal[i]->Add(hist2D);
            hist2D = (NCDataQualityPlot2D*)(tf.Get("DataQuality/"+kDQVars2D[i].name+ncName));
            fDQNearData2DNC[i]->Add(hist2D);
            hist2D = (NCDataQualityPlot2D*)(tf.Get("DataQuality/"+kDQVars2D[i].name+ccName));
            fDQNearData2DCC[i]->Add(hist2D);
          }
        }//end if data

        MSG("NCDataQualityModule", Msg::kInfo) << "  " << fileName
                                                << " -> " << totalName << " "
                                                << ncName << " " << ccName << endl;

        //get each histogram out of the file and add it to the total
        for(int i = 0; i < kDQNumDists; ++i){
          NCDataQualityPlot* hist;

          hist = (NCDataQualityPlot*)(tf.Get("DataQuality/"+kDQVars[i].name+totalName));
          GetDQPlot(i, kDQTotal, datamc, det, kDQNominalIntensity)->Add(hist);
          hist = (NCDataQualityPlot*)(tf.Get("DataQuality/"+kDQVars[i].name+ncName));
          GetDQPlot(i, kDQNC, datamc, det, kDQNominalIntensity)->Add(hist);
          hist = (NCDataQualityPlot*)(tf.Get("DataQuality/"+kDQVars[i].name+ccName));
          GetDQPlot(i, kDQCC, datamc, det, kDQNominalIntensity)->Add(hist);

          if( fileType.Contains("data") || fileType.Contains("blind") ){
            hist = (NCDataQualityPlot*)(tf.Get("LowIntensity/"+kDQVars[i].name+totalName+"LowIntensity"));
            GetDQPlot(i, kDQTotal, datamc, det, kDQLowIntensity)->Add(hist);
            hist = (NCDataQualityPlot*)(tf.Get("LowIntensity/"+kDQVars[i].name+ncName+"LowIntensity"));
            GetDQPlot(i, kDQNC, datamc, det, kDQLowIntensity)->Add(hist);
            hist = (NCDataQualityPlot*)(tf.Get("LowIntensity/"+kDQVars[i].name+ccName+"LowIntensity"));
            GetDQPlot(i, kDQCC, datamc, det, kDQLowIntensity)->Add(hist);
          }

          if(i == kDQEnergy
             && det == Detector::kNear
             && fileName.Contains("data")){
            hist = (NCDataQualityPlot*)(tf.Get("DataQuality/"+kDQVars[i].name+ncName));
            hist->SetDirectory(0);
            name = fileName;
            name.Remove(0, 20);
            name.Remove(7, 17);
            fMonthlySpectraNDNC.push_back(hist);
            fMonthlySpectraNDNC[fMonthlySpectraNDNC.size()-1]->SetName(name+"NC");
            //scale the histogram to 1e19 POT
            const double pot = ((TH1F*)tf.Get("dataPOT"))->Integral();
            fMonthlyNDPOT.push_back(pot);
            MSG("NCDataQualityModule", Msg::kInfo) << "  " << pot << " POT added to " << name << endl;

            hist = (NCDataQualityPlot*)(tf.Get("DataQuality/"+kDQVars[i].name+ccName));
            hist->SetDirectory(0);
            name = fileName;
            name.Remove(0, 20);
            name.Remove(7, 17);
            fMonthlySpectraNDCC.push_back(hist);
          }//end if ND data

        }
      }//end if the right detector and beam type
    }//end loop over files in directory
    delete files;
  }//end loop over beam index vector

  if (datamc==kDQMC && det==Detector::kNear) firstTimeNear=false;
  if (datamc==kDQMC && det==Detector::kFar)  firstTimeFar=false;
}



TString NCDataQualityModule::GetFileTypeString(Detector::Detector_t det,
                                               unsigned int datamc)
{
  // In the case of data and...
  if(datamc == kDQData){
    // Nothing special is happening
    if(!fUseMCAsData && !fUseMockData){
      // In which case use data and blind for near and far respectively
      if(det == Detector::kNear) return "data"; else return "blind";
    }

    // Or we are using a mock data set of some sort in the far detector
    if(fUseMockData && det == Detector::kFar) return fMockDataSet;
  }

  // In all the other cases we use MC files
  return "mc";
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawFarTimingPlots()
{
  //make a canvas for the event locations
  TCanvas* canv3 = new TCanvas("eventSpillTimeCanv", "eventSpillTimeCanv", 150, 10, 900, 600);

  GetDQPlot(kDQDeltaTSpill, kDQTotal, kDQData, Detector::kFar)->Draw();
  GetDQPlot(kDQDeltaTSpill, kDQCC, kDQData, Detector::kFar)->Draw("same");
  GetDQPlot(kDQDeltaTSpill, kDQNC, kDQData, Detector::kFar)->Draw("same");

  TLegend* leg = new TLegend(0.5, 0.5, 0.8, 0.8);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  leg->AddEntry(GetDQPlot(kDQDeltaTSpill, kDQTotal, kDQData, Detector::kFar), "Far Detector Data", "l");
  leg->AddEntry(GetDQPlot(kDQDeltaTSpill, kDQNC, kDQData, Detector::kFar), "NC-like", "l");
  leg->AddEntry(GetDQPlot(kDQDeltaTSpill, kDQCC, kDQData, Detector::kFar), "CC-like", "l");
  leg->Draw();

  MarkPrelim(-20, 50);

  gDirectory->Append(canv3);
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawFarDataQualityPlots()
{
  MSG("NCDataQualityModule", Msg::kInfo) << "far pot - " << fFarPOTTotal
                                         << " in " << fNMCFilesFar << " MC uDSTs"
                                         << endl;


  for(int dist = 0; dist < kDQNumDists; ++dist){
    GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->Write();
    GetDQPlot(dist, kDQCC, kDQData, Detector::kFar)->Write();
    GetDQPlot(dist, kDQNC, kDQData, Detector::kFar)->Write();

//     fDQFarDataTotalLowIntensity[dist]->Write();
//     fDQFarDataCCLowIntensity[dist]->Write();
//     fDQFarDataNCLowIntensity[dist]->Write();

    double scaleFactor=fFarPOTTotal/(fNMCFilesFar*kPotMC);
    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kFar)->Scale(scaleFactor);
    GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar)->Scale(scaleFactor);
    GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar)->Scale(scaleFactor);

    GetDQPlot(dist, kDQTotal, kDQMCOsc, Detector::kFar)->Scale(scaleFactor);
    GetDQPlot(dist, kDQNC, kDQMCOsc, Detector::kFar)->Scale(scaleFactor);
    GetDQPlot(dist, kDQCC, kDQMCOsc, Detector::kFar)->Scale(scaleFactor);

    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kFar)->Write();
//     GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar)->Write();
//     GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar)->Write();

    GetDQPlot(dist, kDQTotal, kDQMCOsc, Detector::kFar)->Write();
//     GetDQPlot(dist, kDQCC, kDQMCOsc, Detector::kFar)->Write();
    GetDQPlot(dist, kDQNC, kDQMCOsc, Detector::kFar)->Write();

    TString canvName = GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->GetName();
    int index = canvName.Index("Data");
    canvName.Remove(index, 4);
    canvName += "Canv";
    TCanvas *canv = new TCanvas(canvName, canvName, 150, 10, 900, 600);
    TLegend *leg = new TLegend(0.75, 0.75, 0.95, 0.95);
    leg->SetBorderSize(0);
    leg->SetFillStyle(0);

    GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->Draw("pe");

    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kFar)->Draw("same ][");
    GetDQPlot(dist, kDQTotal, kDQMCOsc, Detector::kFar)->Draw("same ][");

    leg->AddEntry(GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar), "Far Detector Data", "lp");
    leg->AddEntry(GetDQPlot(dist, kDQTotal, kDQMC, Detector::kFar), "Monte Carlo Total", "l");
    leg->AddEntry(GetDQPlot(dist, kDQTotal, kDQMCOsc, Detector::kFar), "Monte Carlo Total - Oscillated", "l");

    leg->Draw();
    MarkPrelim(GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->GetXaxis()->GetXmax()*0.1,
               GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->GetMaximum()*0.8);

    gDirectory->Append(canv);

    //canvas for NC-like events
    if(fUseNC){
      canvName = GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->GetName();
      index = canvName.Index("Data");
      canvName.Remove(index, 4);
      canvName += "NCCanv";
      TCanvas *canvNC = new TCanvas(canvName, canvName, 150, 10, 900, 600);
      TLegend *leg1 = new TLegend(0.75, 0.75, 0.95, 0.95);
      leg1->SetBorderSize(0);
      leg1->SetFillStyle(0);

      // total strips, r2, pulse height, z
      TString extraDrawOpt("");
      if (dist==kDQTotalStrips  || dist==kDQPulseHeight ||
          dist==kDQRadialSqrVtx)
        extraDrawOpt="][";
      GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar)->Draw(extraDrawOpt);
      GetDQPlot(dist, kDQNC, kDQMCOsc, Detector::kFar)->Draw(extraDrawOpt+"same");
      GetDQPlot(dist, kDQNC, kDQData, Detector::kFar)->Draw("pe same");

      leg1->AddEntry(GetDQPlot(dist, kDQNC, kDQData, Detector::kFar), "Far Detector Data - NC", "p");
      leg1->AddEntry(GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar), "Monte Carlo", "l");
      leg1->AddEntry(GetDQPlot(dist, kDQNC, kDQMCOsc, Detector::kFar), "Monte Carlo - Oscillated ", "l");

      leg1->Draw();
      MarkPrelim(GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar)->GetXaxis()->GetXmax()*0.1,
                 GetDQPlot(dist, kDQNC, kDQMC, Detector::kFar)->GetMaximum()*0.8);

      gDirectory->Append(canvNC);
    }
    if(fUseCC){
      canvName = GetDQPlot(dist, kDQTotal, kDQData, Detector::kFar)->GetName();
      index = canvName.Index("Data");
      canvName.Remove(index, 4);
      canvName += "CCCanv";
      TCanvas *canvCC = new TCanvas(canvName, canvName, 150, 10, 900, 600);
      TLegend *leg2 = new TLegend(0.75, 0.75, 0.95, 0.95);
      leg2->SetBorderSize(0);
      leg2->SetFillStyle(0);

      GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar)->Draw();
      GetDQPlot(dist, kDQCC, kDQData, Detector::kFar)->Draw("pe1same");
      GetDQPlot(dist, kDQCC, kDQMCOsc, Detector::kFar)->Draw("same");
      leg2->AddEntry(GetDQPlot(dist, kDQCC, kDQData, Detector::kFar), "Far Detector Data - CC", "p");
      leg2->AddEntry(GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar), "Monte Carlo", "l");
      leg2->AddEntry(GetDQPlot(dist, kDQCC, kDQMCOsc, Detector::kFar), "Monte Carlo - Oscillated", "l");

      leg2->Draw();
      MarkPrelim(GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar)->GetXaxis()->GetXmax()*0.1,
                 GetDQPlot(dist, kDQCC, kDQMC, Detector::kFar)->GetMaximum()*0.8);

      gDirectory->Append(canvCC);
    }

  }//end loop over dists

  //make plot for executive summary
  TString canvName = "cutsDistsFar";
  TCanvas *canv1 = new TCanvas(canvName, canvName, 150, 10, 900, 600);
  TLegend *leg3 = new TLegend(0.5, 0.6, 0.9, 0.9);
  leg3->SetBorderSize(0);
  leg3->SetFillStyle(0);
  leg3->AddEntry(GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kFar), "Far Detector Data", "lp");
  leg3->AddEntry(GetDQPlot(kDQEventLength, kDQTotal, kDQMCOsc, Detector::kFar), "Oscillated Monte Carlo", "l");

  TPad *pad1 = new TPad("pad1", "", 0.0, 0.5, 0.5, 1.0);
  TPad *pad2 = new TPad("pad2", "", 0.5, 0.5, 1.0, 1.0);
  TPad *pad3 = new TPad("pad3", "", 0.0, 0.0, 0.5, 0.5);
  TPad *pad4 = new TPad("pad4", "", 0.5, 0.0, 1.0, 0.5);

  double maxl = 220.;
  double maxt = 750.;
  double maxe = 50.;

  // the systematic uncertainty bands for the 3.18e20 analysis
  const double elup[50] = {12.663, 12.075,  6.993,  4.237,  3.683,
                           3.620,  3.587,  3.354,  2.932,  2.334,
                           1.873,  1.454,  1.185,  0.989,  0.863,
                           0.761,  0.684,  0.632,  0.581,  0.517,
                           0.473,  0.425,  0.398,  0.271,  0.192,
                           0.218,  0.300,  0.280,  0.252,  0.250,
                           0.217,  0.213,  0.191,  0.182,  0.165,
                           0.149,  0.140,  0.124,  0.111,  0.107,
                           0.095,  0.091,  0.083,  0.077,  0.067,
                           0.061,  0.056,  0.051,  0.015,  0.000};
  const double eldw[50] = {12.663, 12.075,  6.993,  4.237,  3.683,
                           3.620,  3.587,  3.354,  2.932,  2.334,
                           1.873,  1.454,  1.185,  0.989,  0.863,
                           0.761,  0.684,  0.632,  0.581,  0.517,
                           0.473,  0.425,  0.398,  0.271,  0.192,
                           0.218,  0.300,  0.280,  0.252,  0.250,
                           0.217,  0.213,  0.191,  0.182,  0.165,
                           0.149,  0.140,  0.124,  0.111,  0.107,
                           0.095,  0.091,  0.083,  0.077,  0.067,
                           0.061,  0.056,  0.051,  0.015,  0.000};

  const double ntup[50] = {14.531, 51.069,  1.360,  0.001,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000};
  const double ntdw[50] = {14.531, 51.069,  1.360,  0.001,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000};

  const double teup[50] = { 0.001,  0.000,  0.001,  0.001,  0.001,
                            0.001,  0.002,  0.002,  0.002,  0.003,
                            0.003,  0.004,  0.005,  0.010,  0.012,
                            0.019,  0.023,  0.037,  0.062,  0.087,
                            0.131,  0.203,  0.303,  0.455,  0.642,
                            0.903,  1.219,  1.554,  1.806,  1.896,
                            1.707,  1.152,  1.152,  1.030,  0.862,
                            0.439,  0.385,  0.359,  0.342,  0.326,
                            0.317,  0.313,  0.312,  0.299,  0.297,
                            0.303,  0.303,  0.302,  0.290,  0.292};
  const double tedw[50] = { 0.001,  0.000,  0.001,  0.001,  0.001,
                            0.001,  0.002,  0.002,  0.002,  0.003,
                            0.003,  0.004,  0.005,  0.010,  0.012,
                            0.019,  0.023,  0.037,  0.062,  0.087,
                            0.131,  0.203,  0.303,  0.455,  0.642,
                            0.903,  1.219,  1.554,  1.806,  1.896,
                            1.707,  1.152,  1.152,  1.030,  0.862,
                            0.439,  0.385,  0.359,  0.342,  0.326,
                            0.317,  0.313,  0.312,  0.299,  0.297,
                            0.303,  0.303,  0.302,  0.290,  0.292};

  vector<TGraphAsymmErrors *> errors;
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kFar)->GetNbinsX()));
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kFar)->GetNbinsX()));
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kFar)->GetNbinsX()));

  errors[0]->SetName("eventLengthErrGr");
  errors[1]->SetName("numTracksErrGr");
  errors[2]->SetName("trackExtensionErrGr");

  for(int i = 0; i < 3; ++i){
    const double *eup = elup;
    const double *edw = eldw;
    if(i == 1){
      eup = ntup;
      edw = ntdw;
    }
    if(i == 2){
      eup = teup;
      edw = tedw;
    }
    for(int j = 0; j < errors[i]->GetN(); ++j){
      errors[i]->SetPoint(j,
                          GetDQPlot(3+i, kDQTotal, kDQMCOsc, Detector::kFar)->GetBinCenter(j+1),
                          GetDQPlot(3+i, kDQTotal, kDQMCOsc, Detector::kFar)->GetBinContent(j+1));
      errors[i]->SetPointError(j,
                               0.5*GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kFar)->GetBinWidth(j+1),
                               0.5*GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kFar)->GetBinWidth(j+1),
                               edw[j], eup[j]);
    }
  }//end loop to fill error graphs

  TLine *linePlanes = new TLine(60., 0, 60., maxl);
  TLine *lineTracks = new TLine(1., 0, 1., maxt);
  TLine *lineExtension = new TLine(5., 0, 5., maxe);
  TArrow *arrowPlanes = new TArrow(60., 0.75*maxl, 40., 0.75*maxl, 0.02, "|>");
  TArrow *arrowTracks = new TArrow(1., 0.75*maxt, 0.5, 0.75*maxt, 0.02, "|>");
  TArrow *arrowExtension = new TArrow(5., 0.9*maxe, 0., 0.9*maxe, 0.02, "|>");


  linePlanes->SetLineWidth(2);
  lineTracks->SetLineWidth(2);
  lineExtension->SetLineWidth(2);
  linePlanes->SetLineColor(4);
  lineTracks->SetLineColor(4);
  lineExtension->SetLineColor(4);

  arrowPlanes->SetLineWidth(2);
  arrowTracks->SetLineWidth(2);
  arrowExtension->SetLineWidth(2);
  arrowPlanes->SetLineColor(4);
  arrowTracks->SetLineColor(4);
  arrowExtension->SetLineColor(4);
  arrowPlanes->SetFillColor(4);
  arrowTracks->SetFillColor(4);
  arrowExtension->SetFillColor(4);

  GetDQPlot(kDQEventLength, kDQTotal, kDQMCOsc, Detector::kFar)->SetMaximum(maxl);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMCOsc, Detector::kFar)->SetMaximum(maxt);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMCOsc, Detector::kFar)->SetMaximum(maxe);

  pad1->Draw();
  pad2->Draw();
  pad3->Draw();
  pad4->Draw();

  TCanvas* cEvtLength=new TCanvas("event_length_far", "event length far");
  GetDQPlot(kDQEventLength, kDQTotal, kDQMCOsc, Detector::kFar)->GetXaxis()->SetRangeUser(0., 100.);
  GetDQPlot(kDQEventLength, kDQTotal, kDQMCOsc, Detector::kFar)->GetYaxis()->SetRangeUser(0., maxl);
  GetDQPlot(kDQEventLength, kDQTotal, kDQMCOsc, Detector::kFar)->Draw("][");
  errors[0]->Draw("2");
  GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kFar)->Draw("pe1same");
  GetDQPlot(kDQEventLength, kDQCC, kDQMCOsc, Detector::kFar)->Draw("same");
  linePlanes->Draw();
  arrowPlanes->Draw();
  //MarkPrelim(GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kFar)->GetXaxis()->GetXmax()*0.1, maxl*0.8);
  leg3->Draw();
  gDirectory->Append(cEvtLength);

  TCanvas* cNumTracks=new TCanvas("num_tracks_far", "num tracks far");
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMCOsc, Detector::kFar)->GetXaxis()->SetRangeUser(0., 3.);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMCOsc, Detector::kFar)->Draw("][");
  errors[1]->Draw("2");
  GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kFar)->Draw("pe1same");
  GetDQPlot(kDQNumTracks, kDQCC, kDQMCOsc, Detector::kFar)->Draw("same");
  linePlanes->Draw();
  lineTracks->Draw();
  arrowTracks->Draw();
  leg3->Draw();
  //MarkPrelim(GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kFar)->GetXaxis()->GetXmax()*0.1, maxt*0.8);
  gDirectory->Append(cNumTracks);

  //pad3->cd();
  TCanvas* cTrkExt=new TCanvas("trk_ext_far", "track extension far");
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMCOsc, Detector::kFar)->GetXaxis()->SetRangeUser(-12., 20.);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMCOsc, Detector::kFar)->Draw("][");
  errors[2]->Draw("2");
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kFar)->Draw("pe1same");
  GetDQPlot(kDQTrackExtension, kDQCC, kDQMCOsc, Detector::kFar)->Draw("same");
  linePlanes->Draw();
  lineExtension->Draw();
  arrowExtension->Draw();
  leg3->Draw();
  //MarkPrelim(GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kFar)->GetXaxis()->GetXmax()*0.1, maxe*0.8);
  gDirectory->Append(cTrkExt);


  pad4->cd();
  leg3->Draw();

  gDirectory->Append(canv1);

  DrawVertexPlots(Detector::kFar, fFarVtxX, fFarVtxY);
  //DrawVertexPlots(Detector::kFar, fFarVtxXNC, fFarVtxYNC, "NC");
  DrawEventRatePlots(Detector::kFar);
  DrawFarTimingPlots();
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawNearDataQualityPlots()
{
  MSG("NCDataQualityModule", Msg::kInfo) << "near pot - " << fNearPOTTotal
                                         << " in " << fNMCFilesNear << " MC uDSTs"
                                         << endl;

  for(int dist = 0; dist < kDQNumDists; ++dist){

    GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear)->Write();
    GetDQPlot(dist, kDQCC, kDQData, Detector::kNear)->Write();
    GetDQPlot(dist, kDQNC, kDQData, Detector::kNear)->Write();

//     fDQNearDataTotalLowIntensity[dist]->Write();
//     fDQNearDataCCLowIntensity[dist]->Write();
//     fDQNearDataNCLowIntensity[dist]->Write();

    double scaleFactor=fNearPOTTotal/(fNMCFilesNear*kPotMC);
    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kNear)->Scale(scaleFactor);
    GetDQPlot(dist, kDQCC, kDQMC, Detector::kNear)->Scale(scaleFactor);
    GetDQPlot(dist, kDQNC, kDQMC, Detector::kNear)->Scale(scaleFactor);

    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kNear)->Write();
    GetDQPlot(dist, kDQCC, kDQMC, Detector::kNear)->Write();
    GetDQPlot(dist, kDQNC, kDQMC, Detector::kNear)->Write();

    TString canvName = GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear)->GetName();
    int index = canvName.Index("Data");
    canvName.Remove(index, 4);
    canvName += "Canv";
    TCanvas *canv = new TCanvas(canvName, canvName, 150, 10, 900, 600);
    TLegend *leg = new TLegend(0.75, 0.75, 0.95, 0.95);
    leg->SetBorderSize(0);
    leg->SetFillStyle(0);

    GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear)->Draw("pe1");
    GetDQPlot(dist, kDQTotal, kDQMC, Detector::kNear)->Draw("same");
    leg->AddEntry(GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear), "Near Detector Data Total", "lp");
    leg->AddEntry(GetDQPlot(dist, kDQTotal, kDQMC, Detector::kNear), "MC Total", "l");
    if(fUseCC){
      GetDQPlot(dist, kDQCC, kDQData, Detector::kNear)->Draw("pe1same");
      GetDQPlot(dist, kDQCC, kDQMC, Detector::kNear)->Draw("same");
      leg->AddEntry(GetDQPlot(dist, kDQCC, kDQData, Detector::kNear), "Data CC", "p");
      leg->AddEntry(GetDQPlot(dist, kDQCC, kDQMC, Detector::kNear), "MC CC", "l");
    }
    if(fUseNC){
      GetDQPlot(dist, kDQNC, kDQData, Detector::kNear)->Draw("pe1same");
      GetDQPlot(dist, kDQNC, kDQMC, Detector::kNear)->Draw("same");
      leg->AddEntry(GetDQPlot(dist, kDQNC, kDQData, Detector::kNear), "Data NC", "p");
      leg->AddEntry(GetDQPlot(dist, kDQNC, kDQMC, Detector::kNear), "MC NC", "l");
    }

    MarkPrelim(GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear)->GetXaxis()->GetXmax()*0.1,
               GetDQPlot(dist, kDQTotal, kDQData, Detector::kNear)->GetMaximum()*0.8);
    leg->Draw();

    gDirectory->Append(canv);

  }//end loop over dists

  double maxl = 1.25e-5*TMath::Max(GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->GetMaximum(),
                                   GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->GetMaximum());
  double maxt = 1.25e-6*TMath::Max(GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->GetMaximum(),
                                   GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->GetMaximum());
  double maxe = 1.25e-5*TMath::Max(GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->GetMaximum(),
                                   GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->GetMaximum());

  GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->Scale(1.e-5);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->Scale(1.e-6);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->Scale(1.e-5);
  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->Scale(1.e-5);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->Scale(1.e-6);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->Scale(1.e-5);
  GetDQPlot(kDQEventLength, kDQCC, kDQMC, Detector::kNear)->Scale(1.e-5);
  GetDQPlot(kDQNumTracks, kDQCC, kDQMC, Detector::kNear)->Scale(1.e-6);
  GetDQPlot(kDQTrackExtension, kDQCC, kDQMC, Detector::kNear)->Scale(1.e-5);

  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->SetYTitle("Events (#times 10^{5})");
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->SetYTitle("Events (#times 10^{6})");
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->SetYTitle("Events (#times 10^{5})");

  for(int i = 0; i < GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->GetNbinsX(); ++i)
    GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->SetBinError(i+1, 1.e-5*GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->GetBinError(i+1));

  for(int i = 0; i < GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->GetNbinsX(); ++i)
    GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->SetBinError(i+1, 1.e-6*GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->GetBinError(i+1));

  for(int i = 0; i < GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->GetNbinsX(); ++i)
    GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->SetBinError(i+1, 1.e-5*GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->GetBinError(i+1));

  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->SetMaximum(maxl);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->SetMaximum(maxt);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->SetMaximum(maxe);

  //make systematic uncertainty bands

  // PRL error bands
  /*
  double ntup[10] = {40.8293, 191.342, 3.78054, 0.006273,   0,
                     0,       0,       0,       0,          0};
  double ntdw[10] = {37.008,  176.92,  3.57563, 0.00626988, 0,
                     0,       0,       0,       0,          0};
  double elup[50] = {45.5608, 42.911,  22.1133, 16.5154,  15.9571,
                     15.5643, 14.4027, 12.1637,  9.53502,  6.89265,
                      5.043,    3.751,  2.97245, 2.43231,  2.05496,
                      1.81841,  1.6544, 1.60072, 1.50199,  1.6308,
                      3.53188,  3.8482, 3.38007, 2.90701,  2.43652,
                      0.76347,  0,      0,       0.0004,   0,
                     0,0,0,0,0,
                     0,0,0,0,0,
                     0,0,0,0,0,
                     0,0,0,0,0};
  double eldw[50] = {40.4698, 38.8774, 20.4271, 15.2912, 14.9715,
                     14.8539, 13.8931, 11.8686,  9.36871, 6.74724,
                      4.90361, 3.62665, 2.84216, 2.29516, 1.92361,
                      1.6975,  1.53832, 1.48231, 1.39802, 1.51276,
                      3.26844, 3.54244, 3.11487, 2.68013, 2.24106,
                      0.69698, 0,       0,       0.00041, 0,
                     0,0,0,0,0,
                     0,0,0,0,0,
                     0,0,0,0,0,
                     0,0,0,0,0};
  double teup[50] = {0.00181798,0.00155781,0.00125497,0.00188122,0.0013573,
                     0.00268924,0.00327829,0.00455317,0.005519,0.00582066,
                     0.00840255,0.00560191,0.00752972,0.0129454,0.0135563,
                     0.0167458,0.0273413,0.0319429,0.0457396,0.0770695,
                     0.113302,0.181813,0.320604,0.543146,0.931264,
                     1.55251,2.4848,3.77569,5.34883,6.77665,
                     7.69249,6.2814,6.52206,6.14857,5.2589,
                     2.93301,2.48698,2.22001,2.01712,1.88602,
                     1.83302,1.7337,1.69579,1.64394,1.58963,
                     1.5514,1.53539,1.51692,1.55346,1.50479};
  double tedw[50] = {0.00182367,0.00155781,0.00125497,0.00188122,0.0013573,
                     0.00268924,0.00327829,0.00455317,0.005519,0.00582066,
                     0.00840255,0.00560191,0.00752972,0.0129083,0.0135563,
                     0.0167458,0.0273026,0.0319294,0.0457391,0.076818,
                     0.113209,0.1814,0.318877,0.539943,0.920052,
                     1.52943,2.43068,3.65676,5.1071,6.35102,
                     7.07417,5.71773,5.90042,5.54619,4.78795,
                     2.48046,2.11197,1.89213,1.74004,1.63381,
                     1.5828,1.51001,1.48858,1.4472,1.39816,
                     1.37784,1.36899,1.35324,1.381,1.35579};
  */

  const double elup[50] = {51.148,49.002,25.704,19.266,18.943,
                           18.683,17.593,15.076,11.887, 8.496,
                           6.154, 4.612, 3.573, 2.912, 2.454,
                           2.148, 1.954, 1.858, 1.783, 1.921,
                           4.158, 4.478, 3.930, 3.371, 2.794,
                           0.868, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000};
  const double eldw[50] = {51.148, 49.002, 25.704, 19.266, 18.943,
                           18.683, 17.593, 15.076, 11.887,  8.496,
                           6.154,  4.612,  3.573,  2.912,  2.454,
                           2.148,  1.954,  1.858,  1.783,  1.921,
                           4.158,  4.478,  3.930,  3.371,  2.794,
                           0.868,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000};

  const double ntup[50] = {46.800,223.301, 4.484, 0.007, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000,
                           0.000, 0.000, 0.000, 0.000, 0.000};
  const double ntdw[50] = {46.800, 223.301,  4.484,  0.007,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000,
                           0.000,  0.000,  0.000,  0.000,  0.000};

  const double teup[50] = { 0.001, 0.002, 0.002, 0.002, 0.002,
                            0.003, 0.004, 0.004, 0.004, 0.005,
                            0.006, 0.008, 0.012, 0.012, 0.016,
                            0.022, 0.032, 0.043, 0.067, 0.092,
                            0.147, 0.240, 0.415, 0.685, 1.159,
                            1.920, 3.063, 4.589, 6.416, 7.996,
                            8.967, 7.207, 7.513, 6.945, 5.986,
                            3.128, 2.692, 2.403, 2.219, 2.030,
                            1.923, 1.878, 1.846, 1.781, 1.779,
                            1.736, 1.750, 1.731, 1.740, 1.755};
  const double tedw[50] = { 0.001,  0.002,  0.002,  0.002,  0.002,
                            0.003,  0.004,  0.004,  0.004,  0.005,
                            0.006,  0.008,  0.012,  0.012,  0.016,
                            0.022,  0.032,  0.043,  0.067,  0.092,
                            0.147,  0.240,  0.415,  0.685,  1.159,
                            1.920,  3.063,  4.589,  6.416,  7.996,
                            8.967,  7.207,  7.513,  6.945,  5.986,
                            3.128,  2.692,  2.403,  2.219,  2.030,
                            1.923,  1.878,  1.846,  1.781,  1.779,
                            1.736,  1.750,  1.731,  1.740,  1.755};

  vector<TGraphAsymmErrors *> errors;
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->GetNbinsX()));
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->GetNbinsX()));
  errors.push_back(new TGraphAsymmErrors(GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->GetNbinsX()));

  errors[0]->SetName("eventLengthErrGr");
  errors[1]->SetName("numTracksErrGr");
  errors[2]->SetName("trackExtensionErrGr");

  for(int i = 0; i < 3; ++i){
    const double *eup = elup;
    const double *edw = eldw;
    double scale = 1.e-2;
    if(i == 1){
      scale = 1.e-3;
      eup = ntup;
      edw = ntdw;
    }
    if(i == 2){
      eup = teup;
      edw = tedw;
    }
    for(int j = 0; j < errors[i]->GetN(); ++j){
      errors[i]->SetPoint(j,
                          GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kNear)->GetBinCenter(j+1),
                          GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kNear)->GetBinContent(j+1));
      errors[i]->SetPointError(j,
                               0.5*GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kNear)->GetBinWidth(j+1),
                               0.5*GetDQPlot(3+i, kDQTotal, kDQMC, Detector::kNear)->GetBinWidth(j+1),
                               edw[j]*scale, eup[j]*scale);
    }
  }//end loop to fill error graphs

  //make plot for executive summary
  TString canvName = "cutsDists";
  TCanvas *canv1 = new TCanvas(canvName, canvName, 150, 10, 900, 600);
  TLegend *leg1 = new TLegend(0.5, 0.6, 0.9, 0.9);
  leg1->SetBorderSize(0);
  leg1->SetFillStyle(0);
  leg1->AddEntry(GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear), "Near Detector Data", "lp");
  leg1->AddEntry(GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear), "Monte Carlo", "l");
  leg1->AddEntry(GetDQPlot(kDQEventLength, kDQCC, kDQMC, Detector::kNear), "MC CC Background", "f");



  TPad *pad1 = new TPad("pad1", "", 0.0, 0.5, 0.5, 1.0);
  TPad *pad2 = new TPad("pad2", "", 0.5, 0.5, 1.0, 1.0);
  TPad *pad3 = new TPad("pad3", "", 0.0, 0.0, 0.5, 0.5);
  TPad *pad4 = new TPad("pad4", "", 0.5, 0.0, 1.0, 0.5);

  TLine *linePlanes = new TLine(60., 0, 60., maxl);
  TLine *lineTracks = new TLine(1., 0, 1., maxt);
  TLine *lineExtension = new TLine(5., 0, 5., maxe);
  TArrow *arrowPlanes = new TArrow(60., 0.75*maxl, 40., 0.75*maxl, 0.02, "|>");
  TArrow *arrowTracks = new TArrow(1., 0.75*maxt, 0.5, 0.75*maxt, 0.02, "|>");
  TArrow *arrowExtension = new TArrow(5., 0.9*maxe, -5., 0.9*maxe, 0.02, "|>");

  linePlanes->SetLineWidth(2);
  lineTracks->SetLineWidth(2);
  lineExtension->SetLineWidth(2);
  linePlanes->SetLineColor(4);
  lineTracks->SetLineColor(4);
  lineExtension->SetLineColor(4);

  arrowPlanes->SetLineWidth(2);
  arrowTracks->SetLineWidth(2);
  arrowExtension->SetLineWidth(2);
  arrowPlanes->SetLineColor(4);
  arrowTracks->SetLineColor(4);
  arrowExtension->SetLineColor(4);
  arrowPlanes->SetFillColor(4);
  arrowTracks->SetFillColor(4);
  arrowExtension->SetFillColor(4);

  pad1->Draw();
  pad2->Draw();
  pad3->Draw();
  pad4->Draw();

  TCanvas* cEvtLength=new TCanvas("event_length_near", "event length near");
  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->SetRangeUser(0., 100.);
  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->GetYaxis()->SetRangeUser(0., maxl);
  GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->Draw("][");
  errors[0]->Draw("2");
  GetDQPlot(kDQEventLength, kDQCC, kDQMC, Detector::kNear)->Draw("same ][");
  GetDQPlot(kDQEventLength, kDQTotal, kDQData, Detector::kNear)->Draw("pe1same");
  linePlanes->Draw();
  arrowPlanes->Draw();
  //MarkPrelim(GetDQPlot(kDQEventLength, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->GetXmax()*0.1, maxl*0.8);
  leg1->Draw();
  gDirectory->Append(cEvtLength);

  TCanvas* cNumTracks=new TCanvas("num_tracks_near", "num tracks near");
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->SetRangeUser(0., 3.);
  GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->Draw("][");
  errors[1]->Draw("2");
  GetDQPlot(kDQNumTracks, kDQCC, kDQMC, Detector::kNear)->Draw("same ][");
  GetDQPlot(kDQNumTracks, kDQTotal, kDQData, Detector::kNear)->Draw("pe1same");
  lineTracks->Draw();
  arrowTracks->Draw();
  leg1->Draw();
  //MarkPrelim(GetDQPlot(kDQNumTracks, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->GetXmax()*0.1, maxt*0.8);
  gDirectory->Append(cNumTracks);

  //pad3->cd();
  TCanvas* cTrkExt=new TCanvas("trk_ext_near", "track extension near");
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->SetRangeUser(-12., 20.);
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->Draw("][");
  errors[2]->Draw("2");
  GetDQPlot(kDQTrackExtension, kDQCC, kDQMC, Detector::kNear)->Draw("same");
  GetDQPlot(kDQTrackExtension, kDQTotal, kDQData, Detector::kNear)->Draw("pe1same");
  lineExtension->Draw();
  arrowExtension->Draw();
  leg1->Draw();
  //MarkPrelim(GetDQPlot(kDQTrackExtension, kDQTotal, kDQMC, Detector::kNear)->GetXaxis()->GetXmax()*0.1, maxe*0.8);
  gDirectory->Append(cTrkExt);

  //pad4->cd();

  gDirectory->Append(canv1);

  DrawNearEventsPerPOT();
  DrawNearSpectra(fMonthlySpectraNDNC, "NC");
  DrawNearSpectra(fMonthlySpectraNDCC, "CC");
  DrawEventRatePlots(Detector::kNear);
  DrawStabilityPlots(Detector::kNear);
}

//----------------------------------------------------------------------
void NCDataQualityModule::GetNDSectionGraphs(Detector::Detector_t det,
                                             int datamc,
                                             map<int, TGraphErrors*>& spectra,
                                             map<int, TGraphErrors*>& ratios)
{
  NCDataQualityPlot* ptot=GetDQPlot(kDQEnergy, kDQNC, datamc, det);
  for(int i = kDQEnergyQ1; i <= kDQEnergyA2; ++i){
    NCDataQualityPlot* p=GetDQPlot(i, kDQNC, datamc, det);
    spectra[i]=new TGraphErrors(p->GetNbinsX());
    TString name = p->GetName();
    spectra[i]->SetName(name+"Gr");
    ratios[i]=new TGraphErrors(p->GetNbinsX());
    name = p->GetName();
    ratios[i]->SetName(name+"RatioGr");

    const double denom = p->Integral();
    double scale;
    if(denom == 0){
      MSG("NCDataQualityModule", Msg::kInfo) << "PLOT WITH NO INTEGRAL, not scaling" << endl;
      scale = 1;
    } else {
      scale = ptot->Integral()/denom;
    }
    for(int b = 0; b < spectra[i]->GetN(); ++b){
      double x = p->GetBinCenter(b+1);
      double y = p->GetBinContent(b+1);
      double y1 = ptot->GetBinContent(b+1);
      if(i == kDQEnergyQ1) x -= 0.20;
      if(i == kDQEnergyQ2) x -= 0.10;
      if(i == kDQEnergyQ3) x += 0.10;
      if(i == kDQEnergyQ4) x += 0.20;

      if(i == kDQEnergyZ1) x -= 0.25;
      if(i == kDQEnergyZ3) x += 0.25;
      if(i == kDQEnergyA1) x -= 0.25;
      if(i == kDQEnergyA2) x += 0.25;

      spectra[i]->SetPoint(b, x, 1.e-5*scale*y);
      spectra[i]->SetPointError(b, 0., 1.e-5*scale*TMath::Sqrt(y));
      if(y1 > 0) ratios[i]->SetPoint(b, x, scale*y/y1);
      if(y1 > 0 && y > 0) ratios[i]->SetPointError(b, 0., TMath::Sqrt((y*scale/(y1*y1)) + (y1/(y*y*scale*scale))));
    }

  }//end loop to make graphs
}

//----------------------------------------------------------------------
void NCDataQualityModule::FillDoubleRatios(map<int, TGraphErrors*>& nums,
                                           map<int, TGraphErrors*>& denoms,
                                           map<int, TGraphErrors*>& dblratios)
{
  for (map<int, TGraphErrors*>::iterator it=nums.begin();
       it!=nums.end(); ++it) {

    assert(denoms.find(it->first)!=denoms.end());
    TGraphErrors* num=it->second;
    TGraphErrors* denom=denoms.find(it->first)->second;
    //cout << "num=" << num << " denom=" << denom << endl;
    assert(num->GetN()==denom->GetN());
    TGraphErrors* dblratio=new TGraphErrors(num->GetN());
    dblratio->SetName((TString(num->GetName())+"DblRatio").Data());
    dblratios[it->first]=dblratio;
    //cout << "dblratios[" << it->first << "]=" << dblratio << endl;;

    for (int i=0; i<num->GetN(); ++i) {
      double xnum, ynum, xdenom, ydenom;
      num->GetPoint(i, xnum, ynum);
      denom->GetPoint(i, xdenom, ydenom);
      //cout << "ynum=" << ynum << " ydenom=" << ydenom << endl;
      if (ydenom==0) continue;
      double yratio=ynum/ydenom;
      dblratio->SetPoint(i, xnum, yratio);
      double numerr2=SQR(num->GetErrorY(i));
      double denomerr2=SQR(denom->GetErrorY(i));
      dblratio->SetPointError(i, 0, TMath::Sqrt(numerr2 + yratio*denomerr2)/ydenom);
    }
  }
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawStabilityPlot(Detector::Detector_t det,
                                            TString longName,
                                            map<int, TGraphErrors*>& spectra,
                                            map<int, TGraphErrors*>& spectraMC,
                                            map<int, TGraphErrors*>& dblratios,
                                            int minplot, int maxplot,
                                            map<int, TString>& names)
{
  TH1D* forceAxis = new TH1D("ratioForceAxis", ";E_{reco} (GeV); Double Ratio",
                             16, 0., 16);
  forceAxis->SetDirectory(0);
  forceAxis->GetYaxis()->CenterTitle();
  forceAxis->GetYaxis()->SetDecimals();
  // The number 10 below was arrived at by several minutes of trials
  // and error, as the format used by SetNdivisions is maybe the
  // stupidest and most hideous misdesign ever to blight the already
  // disgusting face of ROOT. Can you guess what it does?
  //
  // No, you can't, because it is both unintuitive and impossible to
  // work out from the scant documentation available.
  //
  // (Answer: Don't draw tiny tick marks on the axis.)
  forceAxis->GetYaxis()->SetNdivisions(10);
  forceAxis->GetXaxis()->CenterTitle();

  TCanvas *canv = new TCanvas((longName+"Stability").Data(), "title", 700, 700);
  //TCanvas *canv = new TCanvas("foo");
  TPad *p1 = new TPad("p1", "", 0.0, 0.5, 1.0, 1.0);
  TPad *p2 = new TPad("p2", "", 0.0, 0.0, 1.0, 0.5);

  p1->Draw();
  p2->Draw();

  TLegend *leg = new TLegend(0.62, 0.46, 0.88, 0.9);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);

  p1->cd();
  GetDQPlot(kDQEnergy, kDQNC, kDQData, det)->Draw("hist");
  GetDQPlot(kDQEnergy, kDQNC, kDQMC, det)->Draw("hist same");

  p2->cd();
  forceAxis->Draw();

  for (int v=minplot; v<=maxplot; ++v) {
    p1->cd();
    spectra[v]->Draw("p");
    leg->AddEntry(spectra[v], names[v], "p");
    spectraMC[v]->Draw("p");

    p2->cd();
    //cout << "v=" << v << endl;
    assert(dblratios.find(v)!=dblratios.end());
    //cout << dblratios[v] << endl;
    dblratios[v]->Draw("pe1");
  }

  leg->AddEntry(GetDQPlot(kDQEnergy, kDQNC, kDQData, det),
                "Data Average", "l");
  leg->AddEntry(GetDQPlot(kDQEnergy, kDQNC, kDQMC, det),
                "MC Average", "l");
  p1->cd();
  leg->Draw();

  p2->SetGridy();

  gDirectory->Append(canv);
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawStabilityPlots(Detector::Detector_t det)
{
  //these plots will have 2 pads each, one to overlay the spectra, the
  //other to have the ratios of each section to the total

  //make vectors of graphs to display the data more clearly
  map<int, TGraphErrors*> spectra;
  map<int, TGraphErrors*> ratios;
  // And the same for MC
  map<int, TGraphErrors*> spectraMC;
  map<int, TGraphErrors*> ratiosMC;
  // And the double ratios, because they're the relevant bits
  map<int, TGraphErrors*> dblratios;

  GetNDSectionGraphs(det, kDQData, spectra, ratios);
  GetNDSectionGraphs(det, kDQMC, spectraMC, ratiosMC);

  FillDoubleRatios(ratios, ratiosMC, dblratios);

  for (int i=kDQData; i<=kDQMC; ++i) {
    NCDataQualityPlot* p=GetDQPlot(kDQEnergy, kDQNC, i, det);
    p->Scale(1.e-5);
    p->SetYTitle("Events #times 10^{5}");
    p->GetXaxis()->SetRangeUser(0.,15.);
    p->GetYaxis()->SetDecimals();
  }

  map<int, TString> names;
  names[kDQEnergyQ1]="0 < #theta < #pi/2";
  names[kDQEnergyQ2]="#pi/2 < #theta < #pi";
  names[kDQEnergyQ3]="#pi < #theta < 3#pi/2";
  names[kDQEnergyQ4]="3#pi/2 < #theta < 2#pi";

  names[kDQEnergyZ1]="1.73 < z < 2.73 m";
  names[kDQEnergyZ2]="2.73 < z < 3.73 m";
  names[kDQEnergyZ3]="3.73 < z < 4.73 m";

  names[kDQEnergyA1]="r_{vtx} < 0.5 m";
  names[kDQEnergyA2]="r_{vtx} > 0.5 m";

  DrawStabilityPlot(det, "quadrant", spectra, spectraMC, dblratios,
                    kDQEnergyQ1, kDQEnergyQ4, names);
  DrawStabilityPlot(det, "annuli", spectra, spectraMC, dblratios,
                    kDQEnergyA1, kDQEnergyA2, names);
  DrawStabilityPlot(det, "z", spectra, spectraMC, dblratios,
                    kDQEnergyZ1, kDQEnergyZ3, names);

}


//----------------------------------------------------------------------
void NCDataQualityModule::DrawEventRatePlots(Detector::Detector_t det)
{
  TString detstr = Detector::AsString(det);

  //make event rate vs time plot
  vector<double> date;
  vector<double> dateCC;
  vector<double> dateNC;
  vector<double> pot;
  vector<double> tot;
  vector<double> nc;
  vector<double> cc;
  vector<double> rateTot;
  vector<double> rateNC;
  vector<double> rateCC;

  NCDataQualityPlot* potHist = GetDQPlot(kDQPOTVsTime, kDQTotal, kDQData, det);

  TH1F *forceAxis = new TH1F("forceAxis", ";;Events/10^{18} POT",
                             potHist->GetNbinsX(),
                             potHist->GetXaxis()->GetXmin(),
                             potHist->GetXaxis()->GetXmax());
  forceAxis->GetYaxis()->CenterTitle();

  TH1F *sumPOT = new TH1F("sumPOT", "",
                          potHist->GetNbinsX(),
                          potHist->GetXaxis()->GetXmin(),
                          potHist->GetXaxis()->GetXmax());
  TH1F *sumTot = new TH1F("sumTot", "",
                          potHist->GetNbinsX(),
                          potHist->GetXaxis()->GetXmin(),
                          potHist->GetXaxis()->GetXmax());
  TH1F *sumNC = new TH1F("sumNC", "",
                          potHist->GetNbinsX(),
                          potHist->GetXaxis()->GetXmin(),
                          potHist->GetXaxis()->GetXmax());
  TH1F *sumCC = new TH1F("sumCC", "",
                         potHist->GetNbinsX(),
                         potHist->GetXaxis()->GetXmin(),
                         potHist->GetXaxis()->GetXmax());
  /*
  TH1F *shutDown = new TH1F("shutDown", "",
                            potHist->GetNbinsX(),
                            potHist->GetXaxis()->GetXmin(),
                            potHist->GetXaxis()->GetXmax());
  TH1F *shutDownRate = new TH1F("shutDownRate", "",
                                potHist->GetNbinsX(),
                                potHist->GetXaxis()->GetXmin(),
                                potHist->GetXaxis()->GetXmax());
  */

  for(int i = 0; i < potHist->GetNbinsX(); ++i){
    sumPOT->GetXaxis()->SetBinLabel(i+1, potHist->GetXaxis()->GetBinLabel(i+1));
    sumTot->GetXaxis()->SetBinLabel(i+1, potHist->GetXaxis()->GetBinLabel(i+1));
    sumNC->GetXaxis()->SetBinLabel(i+1, potHist->GetXaxis()->GetBinLabel(i+1));
    sumCC->GetXaxis()->SetBinLabel(i+1, potHist->GetXaxis()->GetBinLabel(i+1));
    forceAxis->GetXaxis()->SetBinLabel(i+1, potHist->GetXaxis()->GetBinLabel(i+1));

    if(potHist->GetBinContent(i+1) > 0.){
      date.push_back(potHist->GetBinCenter(i+1));
      dateCC.push_back(0.1+potHist->GetBinCenter(i+1));
      dateNC.push_back(-0.1+potHist->GetBinCenter(i+1));
      pot.push_back(1.e-6*potHist->GetBinContent(i+1));
      tot.push_back(GetDQPlot(kDQEventsVsTime, kDQTotal, kDQData, det)->GetBinContent(i+1));
      nc.push_back(GetDQPlot(kDQEventsVsTime, kDQNC, kDQData, det)->GetBinContent(i+1));
      cc.push_back(GetDQPlot(kDQEventsVsTime, kDQCC, kDQData, det)->GetBinContent(i+1));

      //fix June 2005 count by hand
      if(i == 1 && det == Detector::kFar) pot[pot.size()-1] += 2.56;

      rateTot.push_back(tot[tot.size()-1]/pot[pot.size()-1]);
      rateNC.push_back(nc[nc.size()-1]/pot[pot.size()-1]);
      rateCC.push_back(cc[cc.size()-1]/pot[pot.size()-1]);

      for(int j = i; j < potHist->GetNbinsX(); ++j){
        sumPOT->Fill(j+0.5, pot[pot.size()-1]);
        sumTot->Fill(j+0.5, tot[pot.size()-1]);
        sumNC->Fill(j+0.5, nc[pot.size()-1]);
        sumCC->Fill(j+0.5, cc[pot.size()-1]);
      }

    }//end if pot in the month
    else{
      sumPOT->SetBinContent(i+1, sumPOT->GetBinContent(i));
      sumTot->SetBinContent(i+1, sumTot->GetBinContent(i));
      sumNC->SetBinContent(i+1, sumNC->GetBinContent(i));
      sumCC->SetBinContent(i+1, sumCC->GetBinContent(i));
    }
  }

  //make graphs of the rates
  TGraphAsymmErrors *totGr = new TGraphAsymmErrors(rateTot.size(), &date[0], &rateTot[0]);
  TGraphAsymmErrors *ncGr = new TGraphAsymmErrors(rateNC.size(), &dateNC[0], &rateNC[0]);
  TGraphAsymmErrors *ccGr = new TGraphAsymmErrors(rateCC.size(), &dateCC[0], &rateCC[0]);

//   TF1 *mean = new TF1("mean"+detector, "pol0");
//   ncGr->Fit(mean, "Q0");
//   TLine *meanLine = new TLine(date[0], mean->GetParameter(0), date[rateNC.size()-1], mean->GetParameter(0));
//   meanLine->SetLineColor(2);
//   meanLine->SetLineWidth(2);

  totGr->SetName("totGr");
  ncGr->SetName("ncGr");
  ccGr->SetName("ccGr");

  //set the errors for the rates
  for (unsigned int i = 0; i < rateTot.size(); ++i){
    double low, high;

    FindUncertainty(tot[i], low, high);
    if(tot[i] > 0) totGr->SetPointError(i, 0, 0, low/pot[i], high/pot[i]);

    FindUncertainty(nc[i], low, high);
    if(nc[i] > 0) ncGr->SetPointError(i, 0, 0, low/pot[i], high/pot[i]);

    FindUncertainty(cc[i], low, high);
    if(cc[i] > 0) ccGr->SetPointError(i, 0, 0, low/pot[i], high/pot[i]);
  }

  TLatex *latex = new TLatex();
  latex->SetName("ShutDown");
  latex->SetTextSize(0.05);
  latex->SetTextColor(4);

  // Box to draw shaded region for shutdown
  const int startShutdown=10, endShutdown=16;
  TBox* shutdownBox=new TBox;
  shutdownBox->SetLineWidth(0);
  shutdownBox->SetFillColor(15);
  shutdownBox->SetFillStyle(3005);
  // Line for start of shutdown
  TLine* shutdownLine1=new TLine;
  TLine* shutdownLine2=new TLine;
  shutdownLine1->SetLineWidth(2);
  shutdownLine1->SetLineColor(15);
  shutdownLine2->SetLineWidth(2);
  shutdownLine2->SetLineColor(15);

  TCanvas *canv = new TCanvas(detstr+"RateCanv", "farRateCanv");
  forceAxis->SetMaximum(3.);
  forceAxis->GetXaxis()->LabelsOption("v");
  forceAxis->Draw();
  canv->Update();
  cout << "rate plot: ymax=" << gPad->GetUymax() << endl;
  shutdownBox->DrawBox(startShutdown, 0, endShutdown, gPad->GetUymax());
  shutdownLine1->DrawLine(startShutdown, 0, startShutdown, gPad->GetUymax());
  shutdownLine2->DrawLine(endShutdown, 0, endShutdown, gPad->GetUymax());

  //totGr->Draw("pe1same");
  ncGr->Draw("pe1same");
  //meanLine->Draw();
  //ccGr->Draw("pe1same");
  /*
  for(int i = 11; i < 17; ++i) shutDownRate->SetBinContent(i, 3.);
  shutDownRate->Draw("same");
  */

  latex->DrawLatex(11, forceAxis->GetMaximum()*0.5, "Shut Down");

  MarkPrelim(forceAxis->GetXaxis()->GetXmax()*0.1, forceAxis->GetMaximum()*0.8);

  TLegend *leg = new TLegend(0.5, 0.5, 0.8, 0.8);
  leg->SetBorderSize(0);
  leg->SetFillStyle(0);
  //leg->AddEntry(totGr, detector+" Detector Data", "p");
  leg->AddEntry(ncGr, "Far Detector NC Events", "lp");
  //leg->AddEntry(ccGr, "Charged Current Events", "p");
  leg->Draw();

  gDirectory->Append(canv);

  TCanvas *canv1 = new TCanvas(detstr + "SummedEvents", "Events");
  sumNC->GetYaxis()->SetRangeUser(0.1, 400);
  sumNC->GetYaxis()->SetLabelSize(0.);
  sumNC->GetYaxis()->SetTitleSize(0.);
  sumNC->GetYaxis()->SetTickLength(0.);
  sumNC->GetXaxis()->LabelsOption("v");
  // Draw it now to set the coord system for the shutdown objects
  sumNC->Draw();
  canv1->Update();
  /*
  for(int i = 11; i < 17; ++i) shutDown->SetBinContent(i, 400);
  shutDown->Draw("same");
  */
  cout << "summed plot: ymax=" << gPad->GetUymax() << endl;
  shutdownBox->DrawBox(startShutdown, 0, endShutdown, gPad->GetUymax());
  shutdownLine1->DrawLine(startShutdown, 0, startShutdown, gPad->GetUymax());
  shutdownLine2->DrawLine(endShutdown, 0, endShutdown, gPad->GetUymax());

  // Draw it again to make it appear on top of the shutdown objects
  sumNC->Draw("same");
  latex->DrawLatex(11, sumNC->GetMaximum()*0.5, "Shut Down");

  //sumTot->Draw("same");
  //sumCC->Draw("same");

  sumPOT->Scale(1.2);
  sumPOT->Draw("same");
  TGaxis *potAxis = new TGaxis(sumTot->GetXaxis()->GetXmax(), 0.,
                               sumTot->GetXaxis()->GetXmax(), 400,
                               // Yay for magic unexplained numbers
                               0., 3.33, 506, "+L");
  potAxis->SetTitle("POT (#times 10^{20})");
  potAxis->SetLabelSize(0.04);
  potAxis->SetLabelFont(42);
  potAxis->SetTitleFont(42);
  potAxis->SetTitleSize(0.055);
  potAxis->SetTitleOffset(0.65);
  potAxis->CenterTitle();
  potAxis->SetLineWidth(2);
  potAxis->SetLineColor(2);
  potAxis->SetTextColor(2);
  potAxis->SetLabelColor(2);
  potAxis->Draw();

  TGaxis *eventAxis = new TGaxis(sumTot->GetXaxis()->GetXmin(), 0.,
                                 sumTot->GetXaxis()->GetXmin(), 400,
                                 0., 400, 506, "");
  eventAxis->SetTitle("Events");
  eventAxis->SetLabelSize(0.04);
  eventAxis->SetLabelFont(42);
  eventAxis->SetTitleFont(42);
  eventAxis->SetTitleSize(0.055);
  eventAxis->SetTitleOffset(0.65);
  eventAxis->CenterTitle();
  eventAxis->SetLineWidth(2);
  eventAxis->Draw();

  MarkPrelim(sumNC->GetXaxis()->GetXmax()*0.1, sumNC->GetMaximum()*0.8);

  TLegend *leg2 = new TLegend(0.5, 0.5, 0.8, 0.8);
  leg2->SetBorderSize(0);
  leg2->SetFillStyle(0);
  leg2->AddEntry(sumPOT, "Protons on Target", "l");
  //leg2->AddEntry(sumTot, "Total "+detector+" Detector Events", "l");
  leg2->AddEntry(sumNC, detstr+" Detector NC Events", "l");
  //leg2->AddEntry(sumCC, "Charged Current Events", "l");
  leg2->Draw();

  gDirectory->Append(canv1);
}

//----------------------------------------------------------------------
void NCDataQualityModule::DrawVertexPlots(Detector::Detector_t det,
                                          std::vector<double>& x,
                                          std::vector<double>& y)
{
  TString detstr = Detector::AsString(det);

  // We want this to be square, but because ROOT is an unutterably
  // stupid and useless bug-ridden mess, we can't use the constructor
  // that takes a "form" argument, like this:
  //
  // TCanvas *canv3 = new TCanvas(detstr+"VtxCanv", detstr+"VtxCanv", 2);
  //
  // because it just ignores this argument in batch mode. Instead, we have to do this:

  TCanvas *canv3 = new TCanvas(detstr+"VtxCanv",
                               detstr+"VtxCanv",
                               500, 500);

  TGraph *vtxGr = new TGraph(x.size(), &x[0], &y[0]);
  vtxGr->SetTitle("");
  vtxGr->SetName(detstr+"VtxYVsXGr");
  vtxGr->GetHistogram()->SetXTitle("x (m)");
  vtxGr->GetHistogram()->SetYTitle("y (m)");
  vtxGr->GetHistogram()->GetXaxis()->CenterTitle();
  vtxGr->GetHistogram()->GetYaxis()->CenterTitle();
  vtxGr->Write();

  float outlineY[] = {1.656, 4.00, 4.00, 1.656, -1.656,
                      -4.00, -4.00, -1.656, 1.656};
  float outlineX[] = {-4.00, -1.656, 1.656, 4.00, 4.00,
                      1.656, -1.656, -4.00, -4.00};
  TGraph *outline = new TGraph(9, outlineX, outlineY);
  outline->SetName("Detector_Outline");
  outline->SetTitle("Detector Outline");
  outline->GetHistogram()->SetXTitle("x (m)");
  outline->GetHistogram()->SetYTitle("y (m)");
  outline->GetHistogram()->GetXaxis()->CenterTitle();
  outline->GetHistogram()->GetYaxis()->CenterTitle();

  float outline1Y[] = {1.449, 3.50, 3.50, 1.449, -1.449,
                       -3.50, -3.50, -1.449, 1.449};
  float outline1X[] = {-3.50, -1.449, 1.449, 3.50, 3.50,
                       1.449, -1.449, -3.50, -3.50};
  TGraph *outline1 = new TGraph(9, outline1X, outline1Y);
  outline1->SetName("FV_Outline");
  outline1->SetTitle("Fiducial Volume");
  TArc *inner = new TArc(0., 0., 0.45);

  outline->Draw("al");
  outline1->Draw("l");
  outline1->SetLineColor(kRed);
  vtxGr->Draw("psame");
  inner->Draw();

  MarkPrelim(-3.5, 4);

  gDirectory->Append(canv3);
}


//----------------------------------------------------------------------
void NCDataQualityModule::DrawNearSpectra(std::vector<NCDataQualityPlot*>& spectra,
                                            TString nccc)
{
  TCanvas *canv1 = new TCanvas("MonthlySpectra"+nccc, "MonthlySpectra"+nccc, 150, 10, 900, 600);
  TLegend *leg1 = new TLegend(0.5, 0.5, 0.75, 0.75);
  leg1->SetBorderSize(0);
  leg1->SetFillStyle(0);
  TLegend *leg2 = new TLegend(0.5, 0.5, 0.75, 0.75);
  leg2->SetBorderSize(0);
  leg2->SetFillStyle(0);

  TString name = "";
  TString months[17] = {"2005-05", "2005-06", "2005-07", "2005-08", "2005-09",
                        "2005-10", "2005-11", "2005-12", "2006-01", "2006-02",
                        "2006-09", "2006-10", "2006-11", "2006-12", "2007-01",
                        "2007-02", "2007-03"};

  if(spectra.empty()) return;

  //make a histogram of the mean for run I and run II
  TH1D *meanRunI = new TH1D("meanRunI"+nccc, "", spectra[0]->GetNbinsX(),
                            spectra[0]->GetXaxis()->GetXmin(),
                            spectra[0]->GetXaxis()->GetXmax());
  TH1D *meanRunII = new TH1D("meanRunII"+nccc, "", spectra[0]->GetNbinsX(),
                             spectra[0]->GetXaxis()->GetXmin(),
                             spectra[0]->GetXaxis()->GetXmax());

  meanRunI->SetXTitle(spectra[0]->GetXaxis()->GetTitle());
  meanRunI->SetYTitle("Events (#times 10^{3})");
  meanRunI->GetXaxis()->CenterTitle();
  meanRunI->GetYaxis()->CenterTitle();

  meanRunII->SetXTitle(spectra[0]->GetXaxis()->GetTitle());
  meanRunII->SetYTitle("Events (#times 10^{3})");
  meanRunII->GetXaxis()->CenterTitle();
  meanRunII->GetYaxis()->CenterTitle();

  //make a vector of tgraphs for run I and run II
  vector<TGraphErrors *> runI;
  vector<TGraphErrors *> runII;

  for (int m = 0; m < 17; ++m){
    for (unsigned int i = 0; i < spectra.size(); ++i){
      name = spectra[i]->GetName();
      if(name.Contains(months[m])){
        if(m < 10){
          runI.push_back(new TGraphErrors(spectra[i]->GetNbinsX()));
          runI[runI.size()-1]->SetName(months[m]);
          for(int j = 0; j < spectra[i]->GetNbinsX(); ++j){
            runI[runI.size()-1]->SetPoint(j, 0.1*(m-5)+0.05+spectra[i]->GetBinCenter(j+1),
                                          1.e-3*(1.e7/fMonthlyNDPOT[i])*spectra[i]->GetBinContent(j+1));
            runI[runI.size()-1]->SetPointError(j, 0.,
                                               1.e-3*(1.e7/fMonthlyNDPOT[i])*spectra[i]->GetBinError(j+1));
          }
          meanRunI->Add(spectra[i], 1.e-3*(1.e7/fMonthlyNDPOT[i])*0.1);
        }//end if run I
        else{
          runII.push_back(new TGraphErrors(spectra[i]->GetNbinsX()));
          runII[runII.size()-1]->SetName(months[m]);
          for(int j = 0; j < spectra[i]->GetNbinsX(); ++j){
            runII[runII.size()-1]->SetPoint(j, 0.1*(m-13)+0.05+spectra[i]->GetBinCenter(j+1),
                                            1.e-3*(1.e7/fMonthlyNDPOT[i])*spectra[i]->GetBinContent(j+1));
            runII[runII.size()-1]->SetPointError(j, 0.,
                                                 1.e-3*(1.e7/fMonthlyNDPOT[i])*spectra[i]->GetBinError(j+1));
          }
          meanRunII->Add(spectra[i], 1.e-3*(1.e7/fMonthlyNDPOT[i])*0.142857);
        }//end if run II
      }//end if on the right month
    }//end loop over spectra
  }//end loop over months

  //get skzp error envelope
  double skzpErrY[15] = {0.71, 0.43, 0.38, 0.39, 0.34,
                         0.24, 0.14, 0.09, 0.07, 0.06,
                         0., 0., 0., 0., 0.};

  TGraphAsymmErrors *meanRunIGr = new TGraphAsymmErrors(15);
  TGraphAsymmErrors *meanRunIIGr = new TGraphAsymmErrors(15);

  meanRunIGr->SetTitle("");
  meanRunIIGr->SetTitle("");

  meanRunIGr->SetName("meanRunIGr");
  meanRunIIGr->SetName("meanRunIIGr");

  for(int i = 0; i < 15; ++i){
    meanRunIGr->SetPoint(i, meanRunI->GetBinCenter(i+1),
                         meanRunI->GetBinContent(i+1));
    meanRunIGr->SetPointError(i, 0.5, 0.5, skzpErrY[i], skzpErrY[i]);

    meanRunIIGr->SetPoint(i, meanRunII->GetBinCenter(i+1),
                         meanRunII->GetBinContent(i+1));
    meanRunIIGr->SetPointError(i, 0.5, 0.5, skzpErrY[i], skzpErrY[i]);
  }

  meanRunIGr->GetHistogram()->GetXaxis()->SetTitle("E_{vis} (GeV)");
  meanRunIIGr->GetHistogram()->GetXaxis()->SetTitle("E_{vis} (GeV)");

  meanRunIGr->GetHistogram()->GetYaxis()->SetTitle("Events/10^{19} POT");
  meanRunIIGr->GetHistogram()->GetYaxis()->SetTitle("Events/10^{19} POT");

  meanRunIGr->GetHistogram()->GetXaxis()->CenterTitle();
  meanRunIIGr->GetHistogram()->GetXaxis()->CenterTitle();

  meanRunIGr->GetHistogram()->GetYaxis()->CenterTitle();
  meanRunIIGr->GetHistogram()->GetYaxis()->CenterTitle();

  TPad *pad1 = new TPad("pad1", "", 0.0, 0.5, 1.0, 1.0);
  TPad *pad2 = new TPad("pad2", "", 0.0, 0.0, 1.0, 0.5);

  pad1->Draw();
  pad2->Draw();

  pad1->cd();
  meanRunIGr->GetHistogram()->GetYaxis()->SetRangeUser(0.1, 12);
  meanRunIGr->Draw("pa2");
  meanRunI->Draw("same");
  for (unsigned int i = 0; i < runI.size(); ++i){
    runI[i]->Draw("pe1same");
    leg1->AddEntry(runI[i], runI[i]->GetName(), "lp");
  }
  leg1->AddEntry(meanRunI, "ND Run I Mean - "+nccc, "l");
  leg1->SetNColumns(2);
  leg1->Draw();

  pad2->cd();
  meanRunIIGr->GetHistogram()->GetYaxis()->SetRangeUser(0.1, 12);
  meanRunIIGr->Draw("pa2");
  meanRunII->Draw("same");
  for (unsigned int i = 0; i < runII.size(); ++i){
    runII[i]->Draw("pe1same");
    leg2->AddEntry(runII[i], runII[i]->GetName(), "lp");
  }
  leg2->AddEntry(meanRunII, "ND Run II Mean - "+nccc, "l");
  leg2->SetNColumns(2);
  leg2->Draw();

  gDirectory->Append(canv1);
}



//----------------------------------------------------------------------
void NCDataQualityModule::DrawNearEventsPerPOT()
{

  for(int dist = 0; dist < kDQ2DNumDists; ++dist){
    fDQNearData2DTotal[dist]->Write();
    fDQNearData2DCC[dist]->Write();
    fDQNearData2DNC[dist]->Write();
  }

  //make vectors of profile histograms
  vector<TProfile *> profTot;
  vector<TProfile *> profNC;
  vector<TProfile *> profCC;

  //and of TF1
  vector<TF1 *> fitTot;
  vector<TF1 *> fitNC;
  vector<TF1 *> fitCC;
  TString name = "";

  for (unsigned int i = 0; i < fDQNearData2DTotal.size(); ++i){
    profTot.push_back(fDQNearData2DTotal[i]->ProfileX());
    name = fDQNearData2DTotal[i]->GetName();
    fitTot.push_back(new TF1("fit"+name, "pol0"));
    profTot[i]->Fit(fitTot[i], "Q0");
    profNC.push_back(fDQNearData2DNC[i]->ProfileX());
    name = fDQNearData2DNC[i]->GetName();
    fitNC.push_back(new TF1("fit"+name, "pol0"));
    profNC[i]->Fit(fitNC[i], "Q0");
    profCC.push_back(fDQNearData2DCC[i]->ProfileX());
    name = fDQNearData2DCC[i]->GetName();
    fitCC.push_back(new TF1("fit"+name, "pol0"));
    profCC[i]->Fit(fitCC[i], "Q0");
  }

  //draw the profiles
  TCanvas *canv1 = new TCanvas("totalEventsPerPOT", "totalEventsPerPOT", 150, 10, 900, 600);
  TLegend *leg1 = new TLegend(0.75, 0.75, 0.95, 0.95);
  leg1->SetBorderSize(0);
  leg1->SetFillStyle(0);

  profTot[kDQ2DEventsPerPOT]->Draw();
  profTot[kDQ2DEventsPerPOTLT05GeV]->Draw("same");
  profTot[kDQ2DEventsPerPOTLT1GeV]->Draw("same");
  profTot[kDQ2DEventsPerPOTGT1GeV]->Draw("same");
  fitTot[kDQ2DEventsPerPOTLT05GeV]->Draw();
  fitTot[kDQ2DEventsPerPOTLT1GeV]->Draw();
  fitTot[kDQ2DEventsPerPOTGT1GeV]->Draw();
  fitTot[kDQ2DEventsPerPOT]->Draw();

  leg1->AddEntry(profTot[kDQ2DEventsPerPOT], "All Events", "p");
  leg1->AddEntry(profTot[kDQ2DEventsPerPOTLT05GeV], "0 < E_{vis} < 0.5 GeV", "p");
  leg1->AddEntry(profTot[kDQ2DEventsPerPOTLT1GeV], "0.5 < E_{vis} < 1.0 GeV", "p");
  leg1->AddEntry(profTot[kDQ2DEventsPerPOTGT1GeV], "1 GeV < E_{vis}", "p");
  leg1->Draw();

  gDirectory->Append(canv1);

  TCanvas *canv2 = new TCanvas("ncEventsPerPOT", "ncEventsPerPOT", 150, 10, 900, 600);
  TLegend *leg2 = new TLegend(0.75, 0.75, 0.95, 0.95);
  leg2->SetBorderSize(0);
  leg2->SetFillStyle(0);

  profNC[kDQ2DEventsPerPOT]->Draw();
  profNC[kDQ2DEventsPerPOTLT05GeV]->Draw("same");
  profNC[kDQ2DEventsPerPOTLT1GeV]->Draw("same");
  profNC[kDQ2DEventsPerPOTGT1GeV]->Draw("same");
  fitNC[kDQ2DEventsPerPOTLT05GeV]->Draw();
  fitNC[kDQ2DEventsPerPOTLT1GeV]->Draw();
  fitNC[kDQ2DEventsPerPOTGT1GeV]->Draw();
  fitNC[kDQ2DEventsPerPOT]->Draw();

  leg2->AddEntry(profNC[kDQ2DEventsPerPOT], "All Events", "p");
  leg2->AddEntry(profNC[kDQ2DEventsPerPOTLT05GeV], "0 < E_{vis} < 0.5 GeV", "p");
  leg2->AddEntry(profNC[kDQ2DEventsPerPOTLT1GeV], "0.5 < E_{vis} < 1.0 GeV", "p");
  leg2->AddEntry(profNC[kDQ2DEventsPerPOTGT1GeV], "1 GeV < E_{vis}", "p");
  leg2->Draw();

  gDirectory->Append(canv2);

  TCanvas *canv3 = new TCanvas("ccEventsPerPOT", "ccEventsPerPOT", 150, 10, 900, 600);
  TLegend *leg3 = new TLegend(0.75, 0.75, 0.95, 0.95);
  leg3->SetBorderSize(0);
  leg3->SetFillStyle(0);

  profCC[kDQ2DEventsPerPOT]->Draw();
  profCC[kDQ2DEventsPerPOTLT05GeV]->Draw("same");
  profCC[kDQ2DEventsPerPOTLT1GeV]->Draw("same");
  profCC[kDQ2DEventsPerPOTGT1GeV]->Draw("same");
  fitCC[kDQ2DEventsPerPOTLT05GeV]->Draw();
  fitCC[kDQ2DEventsPerPOTLT1GeV]->Draw();
  fitCC[kDQ2DEventsPerPOTGT1GeV]->Draw();
  fitCC[kDQ2DEventsPerPOT]->Draw();

  leg3->AddEntry(profCC[kDQ2DEventsPerPOT], "All Events", "p");
  leg3->AddEntry(profCC[kDQ2DEventsPerPOTLT05GeV], "0 < E_{vis} < 0.5 GeV", "p");
  leg3->AddEntry(profCC[kDQ2DEventsPerPOTLT1GeV], "0.5 < E_{vis} < 1.0 GeV", "p");
  leg3->AddEntry(profCC[kDQ2DEventsPerPOTGT1GeV], "1 GeV < E_{vis}", "p");
  leg3->Draw();

  gDirectory->Append(canv3);
}


void NCDataQualityModule::MarkPrelim(double x, double y)
{
  static TLatex* latex = 0;

  if(!latex){
    latex = new TLatex();
    latex->SetName("Prelim");
    latex->SetTextSize(0.05);
    latex->SetTextColor(kBlue);
   }

  static const TString text = TString::Format("MINOS Preliminary: %3.2lf #times 10^{20} POT",
                                              fFarPOTTotal/1e20);

  latex->DrawLatex(x, y, text);
}


//-----------------------------------------------------------------------
void NCDataQualityModule::FindUncertainty(double value,
                                          double& uncertaintyLow,
                                          double& uncertaintyHigh)
{
  const double poissonErrorsUp[] = {1.29,  2.75,  4.25,  5.30,  6.78,
                                    7.81,  9.28, 10.30, 11.32, 12.79,
                                   13.81, 14.82, 16.29, 17.30, 18.32,
                                   19.32, 20.80, 21.81, 22.82, 23.82};
  const double poissonErrorsDown[] = {0.00,  0.37,  0.74,  1.10,  2.34,
                                      2.75,  3.82,  4.25,  5.30,  6.33,
                                      6.78,  7.81,  8.83,  9.28, 10.30,
                                     11.32, 12.33, 12.79, 13.81, 14.82};

  if(value > 19){
    uncertaintyLow = TMath::Sqrt(value);
    uncertaintyHigh = TMath::Sqrt(value);
  }
  else{
    uncertaintyHigh = poissonErrorsUp[int(value)]-value;
    uncertaintyLow = value - poissonErrorsDown[int(value)];
  }
}

---