{
#include <TH1D.h>
#include <TFile.h>
#include <TUnixSystem.h>
#include <TGraphErrors.h>
#include <TGraph.h>
#include <TF1.h>

char* inputName = "/home5/dwa123/offline/work/workspace/V01-03-03/output/str30/cut15/merged/DistanceMergedSingle_LogGauss_Split.root";
TFile* inputFile = new TFile(inputName);

char* outputName = "/home5/dwa123/offline/work/workspace/V01-03-03/output/str30/cut15/merged/geometryLogGaussAlt_Split.root";
TFile* outputFile = new TFile(outputName, "RECREATE");

double doms[59];
double zenMeans[59][4];
double errDoms[59];
double zenErrMeans[59][4];
double zenRmsMeans[59][4];
double distMeans[59][4];
double distErrMeans[59][4];
double distRmsMeans[59][4];

cout << "Began looping over DOMs." << endl;
for (int iDOM = 0; iDOM < 59; iDOM++)  {

  for (int iSplit = 0; iSplit < 4; iSplit++)  {

    //construct the name for the histograms associated with this particular
    //split group
    char zenName[255];
    char distName[255];

    double dSplit = (double) iSplit;
    sprintf(zenName, "DOM %d,%d, zen%g", (iDOM+1), (iDOM+2), 1.0-dSplit*0.05);
    sprintf(distName, "DOM %d,%d, dist%d", iDOM+1, iDOM+2, iSplit*10);

    //extract the histograms from the file
    TH1D* zenHist = (TH1D*) inputFile->Get(zenName);
    TH1D* distHist = (TH1D*) inputFile->Get(distName);

    if (zenHist)  {

      TF1* fitFun = zenHist->GetFunction("logN");
      double rms = zenHist->GetRMS();

      double sigma = fitFun->GetParameter(0);
      double errSigma = fitFun->GetParError(0);
      double mu = fitFun->GetParameter(2);
      double errMu = fitFun->GetParError(2);
      double offset = fitFun->GetParameter(1);
      double errOffset = fitFun->GetParError(1);
      double logMean = exp(mu - sigma*sigma);
      double errLogMean = logMean * sqrt(errMu*errMu + 4*sigma*sigma*errSigma*errSigma);
      double logFit = logMean - offset;
      double errLogFit = TMath::Abs(logFit) * sqrt(errLogMean*errLogMean/(logMean*logMean) + errOffset*errOffset/(offset*offset));

      zenMeans[iDOM][iSplit] = logFit;
      zenErrMeans[iDOM][iSplit] = errLogFit;
      doms[iDOM] = iDOM + 1;
      errDoms[iDOM] = 0;
      zenRmsMeans[iDOM][iSplit] = rms;
    }else {

      cerr << "Failed to recover zenith histogram #" << iDOM << ", " << iSplit
           << endl;

      doms[iDOM] = iDOM + 1;
      zenMeans[iDOM][iSplit] = 0;
      errDoms[iDOM] = 0;
      zenErrMeans[iDOM][iSplit] = 0;
    }

    if (distHist)  {

      TF1* fitFun = distHist->GetFunction("logN");
      double rms = distHist->GetRMS();

      double sigma = fitFun->GetParameter(0);
      double errSigma = fitFun->GetParError(0);
      double mu = fitFun->GetParameter(2);
      double errMu = fitFun->GetParError(2);
      double offset = fitFun->GetParameter(1);
      double errOffset = fitFun->GetParError(1);
      double logMean = exp(mu - sigma*sigma);
      double errLogMean = logMean * sqrt(errMu*errMu + 4*sigma*sigma*errSigma*errSigma);
      double logFit = logMean - offset;
      double errLogFit = TMath::Abs(logFit) * sqrt(errLogMean*errLogMean/(logMean*logMean) + errOffset*errOffset/(offset*offset));

      distMeans[iDOM][iSplit] = logFit;
      distErrMeans[iDOM][iSplit] = errLogFit;
      distRmsMeans[iDOM][iSplit] = rms;
    }else {

      cerr << "Failed to recover distance histogram #" << iDOM << ", " << iSplit
           << endl;

      doms[iDOM] = iDOM + 1;
      distMeans[iDOM][iSplit] = 0;
      errDoms[iDOM] = 0;
      distErrMeans[iDOM][iSplit] = 0;
    }
  }//end loop over Split groups
}

double total = 0.0;
for (int iDOM = 0; iDOM < 58; iDOM++)  {

  total += zenMeans[iDOM][0];
}
double mean = total /59;
cout << "mean = " << mean << endl;

//loop over all the split groups, and create a geometry summary for each
for (int iSplit = 0; iSplit < 4; iSplit++)  {

//  cout << "Split #" << iSplit << endl;

  //create a double to hold the value of the split index
  double dSplit = (double) iSplit;

  //create fresh arrays to hold the actual plotting data
  double zenPlot[59];
  double distPlot[59];
  double zenPlotErr[59];
  double distPlotErr[59];
  double zenPlotRMS[59];
  double distPlotRMS[59]

  //fill the graphing arrays
  for (int iDOM = 0; iDOM < 59; iDOM++)  {

    zenPlot[iDOM] = zenMeans[iDOM][iSplit];
    zenPlotErr[iDOM] = zenErrMeans[iDOM][iSplit];
    zenPlotRMS[iDOM] = zenRmsMeans[iDOM][iSplit];
    distPlot[iDOM] = distMeans[iDOM][iSplit];
    distPlotErr[iDOM] = distErrMeans[iDOM][iSplit];
    distPlotRMS[iDOM] = distRmsMeans[iDOM][iSplit];
  }

  //initialize and fill the names of the four graphs about to be created
  char zenErrName[255];
  char distErrName[255];
  char zenRmsName[255];
  char distRmsName[255];
  sprintf(zenErrName, "summaryErr_zen%g", 1.0-dSplit*0.05);
  sprintf(distErrName, "summaryErr_dist%d", 10*iSplit);
  sprintf(zenRmsName, "summaryRMS_zen%g", 1.0-dSplit*0.05);
  sprintf(distRmsName, "summaryRMS_dist%g", 10*iSplit);

  //write zenith splitting geometries
  TGraphErrors* summaryErrZen = new TGraphErrors(59, doms, zenPlot, errDoms, 
                                                 zenPlotErr);
  TF1 func("p0", "[0]", 1, 30);
  func.SetParameter(0, 0);

  summaryErr->Fit("p0", "RM");

  summaryErrZen->SetName(zenErrName);
  summaryErrZen->SetTitle("Differences in Distances Between Adjacent Pairs(m)");
  summaryErrZen->GetXaxis()->SetTitle("Pair Number");
  summaryErrZen->GetYaxis()->SetTitle("Difference in Distances(m)");

  summaryErrZen->Write();

  TGraphErrors* summaryRMSZen = new TGraphErrors(59, doms, zenPlot, errDoms, 
                                                 zenPlotRMS);
  TF1 func1("p1", "[0]", 1, 30);
  func1.SetParameter(0, 0);

  summaryRMSZen->Fit("p1", "RM");
  summaryRMSZen->SetName(zenRmsName);
  summaryRMSZen->SetTitle("Difference Between Expected and Measured Distances");
  summaryRMSZen->GetXaxis()->SetTitle("Pair Number");
  summaryRMSZen->GetYaxis()->SetTitle("Difference in Distances (m)");

  summaryRMSZen->Write();

  //write distance splitting geometries
  TGraphErrors* summaryErrDist = new TGraphErrors(59, doms, distPlot, errDoms,
                                                  distPlotErr);
  TF1 func2("p2", "[0]", 1, 30);
  func1.SetParameter(0, 0);

  summaryErrDist->Fit("p2", "RM");
  summaryErrDist->SetName(distErrName);
  summaryErrDist->SetTitle("Differences in Distances Between Adjacent Pairs(m)");
  summaryErrDist->GetXaxis()->SetTitle("Pair Number");
  summaryErrDist->GetYaxis()->SetTitle("Difference in Distances(m)");

  summaryErrDist->Write();

  TGraphErrors* summaryRMSDist = new TGraphErrors(59, doms, distPlot, errDoms, 
                                                 distPlotRMS);
  TF1 func3("p3", "[0]", 1, 30);
  func3.SetParameter(0, 0);

  summaryRMSDist->Fit("p3", "RM");
  summaryRMSDist->SetName(zenRmsName);
  summaryRMSDist->SetTitle("Difference Between Expected and Measured Distances");
  summaryRMSDist->GetXaxis()->SetTitle("Pair Number");
  summaryRMSDist->GetYaxis()->SetTitle("Difference in Distances (m)");
}

outputFile->Close(); 
}