// This script is meant to be "standard" analysis code for the July-August 2012 beam test at TRIUMF.
// You call it using ROOT as follows, from the bash prompt:
// > root -l -b -q 'standard_analyzed.C+("path_to_dir")'
// Where path_to_dir is where are located the C1, C2 and C3.root files for the run that you wish to analyse.
// The output is a .root file in the same directory containing some sample graphs and a tree with
// the processed variables.  The processed variables are described in the section which creates tree branches.

#define N_TRACEPLOTS 0 // Number of events for which to produce sample graphs.
#define DEBUG 0 // Set to 1 to see debug messages.
#define N_SAMPLES 20002 // The number of voltage samples in each oscilloscope trace.

#include "TVectorD.h"
#include "TFile.h"
#include "TTree.h"
#include <iostream>
#include <vector>
#include <algorithm>
#include "TH1F.h"
#include "TCanvas.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TString.h"
#include "TGraph.h"
#include "TDirectory.h"
#include "TPad.h"
#include "TMath.h"

using namespace std;

// These are function prototypes, the actual bodies of the functions are below the main standard_analysis function.
vector<Int_t> clu_thresh_over_average(TVectorD ampl_tv, Int_t n_frames, Float_t thresh);
Float_t get_scintilator_distance(TString path_to_dir);
TVectorD smoothed_tv(TVectorD ampl_tv, Int_t n_frames);

// This is the main function which is run.
void standard_analysis(TString path_to_dir)
{
  // Open files
  TFile f1(path_to_dir+"/C1.root");
  TFile f2(path_to_dir+"/C2.root");
  TFile f3(path_to_dir+"/C3.root");

  TFile outfile(path_to_dir+"/standard_analysed.root","RECREATE");
 
  // Get trees and set addresses for the relevant branches.
  f1.cd();
  TTree* C1 = (TTree*)gDirectory->Get("C1");
  TVectorD *C1_tv=0;
  C1->SetBranchAddress("ampl_tv",&C1_tv);

  f2.cd();
  TTree* C2 = (TTree*)gDirectory->Get("C2");
  TVectorD *C2_tv=0;
  C2->SetBranchAddress("ampl_tv",&C2_tv);

  f3.cd();
  TTree* C3 = (TTree*)gDirectory->Get("C3");
  TVectorD *C3_tv=0;
  C3->SetBranchAddress("ampl_tv",&C3_tv);

  if(DEBUG == 1)
    {
      cout << "C1 entries: " << C1->GetEntries() << endl;
      cout << "C2 entries: " << C2->GetEntries() << endl;
      cout << "C3 entries: " << C3->GetEntries() << endl;
    }

  // Before looping, create output tree and branches.
  outfile.cd();
  TTree *outtree = new TTree("standard_analyzed","Tree with processed quantities.");

  // These are variables that are used during various calculations.
  Float_t scintillator_distance = get_scintilator_distance(path_to_dir);
  const Int_t n_C1_pedestalframes = 1000; // For TOF pedestal determination
  Float_t TOF_threshold = -0.4; // For TOF calculation

  // These variables store the baselines and RMS of the empty events, used to set thresholds for the charge integrations.
  Float_t prevIntegralC2 = 0,newIntegralC2 = 0; // For DiffUncorr
  Float_t prevIntegralC3 = 0,newIntegralC3 = 0; // For DiffUncorr
  Double_t baseline2 = 0, baseline3 = 0; // In volts, baseline of smoothed empty signal
  Double_t prevRMS2 = 0, prevRMS3 = 0; // In volts, RMS of smoothed empty signals
  Float_t RMSMult = -5.0; // Fluctuation from baseline in units of RMS to trigger signal threshold.
  Int_t StartIntegral2 = -99, StartIntegral3 = -99; // Where the signal threshold was crossed.
  Int_t prespace = 100; // Frames to skip back from signal threshold for charge integral.
  // Number of frames over which to smooth for signal threshold crossing, 100 samples means I am smoothing over 5 nanoseconds.
  const Int_t n_frames_smooth = 100; 
  TVectorD smoothed2(N_SAMPLES),smoothed3(N_SAMPLES);
  TVectorD sub(n_C1_pedestalframes),sub2,sub3; // Subvectors for the TOF and charge integrals.
  const Int_t integral_duration = 650; // 13000 samples is 650 ns
  Int_t endintegral2 = 0, endintegral3 = 0;
  Double_t sub_baseline2 = 0, sub_baseline3 = 0;
  Double_t Charge2 = 0, Charge3 = 0; // Not scaled to pC, this is just sums of sample values.

  Float_t input_impedance = 75.0; // Ohms
  Float_t sample_spacing = 50.0; // picoseconds

  Float_t cc_threshold = 0.005; // For cluster counting algorithm (in volts)
  Int_t cc_n_frames = 250; // For cluster counting algorithm (in frames)

  // These are variables which will be linked to branches in outtree.  Make sure
  // that these are set to valid or "nonsense" values in every iteration over an
  // event, because *all* of them are written in each entry in the tree.

  Int_t rawTOF = 0,nTOF = 0;
  Double_t TOF = 0,Beta = 0;
  Bool_t SignalPresent2 = kFALSE, SignalPresent3 = kFALSE; // Flags for threshold algorithm
  Double_t Charge2_pC = 0,Charge3_pC = 0,Pedestal1 = 0,Pedestal2_pC = 0,Pedestal3_pC = 0;
  Double_t RawCharge2_pC = 0,RawCharge3_pC = 0,MinimumVal2 = 0,MinimumVal3 = 0;
  Double_t FracCharge2 = 0, FracCharge3 = 0;
  Double_t DiffUncorr2 = 0,DiffUncorr3 = 0;

  vector<Int_t> Clusters2;
  vector<Int_t> Clusters3;

  if(DEBUG == 1){cout << "Creating branches on outtree." << endl;}

  outtree->Branch("SignalPresent2",&SignalPresent2,"SignalPresent2/O"); // If channel 2 had a threshold crossing
  outtree->Branch("SignalPresent3",&SignalPresent3,"SignalPresent3/O"); // If channel 3 had a threshold crossing
  outtree->Branch("rawTOF",&rawTOF,"rawTOF/I"); // Unscaled TOF calculation (i.e. in TDC counts)
  outtree->Branch("TOF",&TOF,"TOF/D"); // Scaled TOF calculation (in nanoseconds)
  outtree->Branch("nTOF",&nTOF,"nTOF/I"); // Number of TOF thresholds found, should be 4 for good triggers, 0 for Uncorr
  outtree->Branch("Beta",&Beta,"Beta/D"); // Speed of particle calculated from TOF (distance needs to be adjusted per run)
  outtree->Branch("Charge2_pC",&Charge2_pC,"Charge2_pC/D"); // Integrated charge of 20um chamber signal (in picoCoulombs)
  outtree->Branch("Charge3_pC",&Charge3_pC,"Charge3_pC/D"); // Integrated charge of 25um chamber signal (in picoCoulombs)
  outtree->Branch("Pedestal1",&Pedestal1,"Pedestal1/D"); // Average of first 100 samples of trigger signal
  outtree->Branch("Pedestal2_pC",&Pedestal2_pC,"Pedestal2_pC/D"); // Integrated signal of uncorrelated signals (in pC)
  outtree->Branch("Pedestal3_pC",&Pedestal3_pC,"Pedestal3_pC/D"); // Integrated signal of uncorrelated signals (in pC)
  outtree->Branch("RawCharge2_pC",&RawCharge2_pC,"RawCharge2_pC/D"); // Same as Charge2 but without pedestal subtraction (in picoCoulombs)
  outtree->Branch("RawCharge3_pC",&RawCharge3_pC,"RawCharge3_pC/D"); // Same as Charge3 but without pedestal subtraction (in picoCoulombs)
  outtree->Branch("MinimumVal2",&MinimumVal2,"MinimumVal2/D"); // Smallest voltage sample of 20um chamber signal
  outtree->Branch("MinimumVal3",&MinimumVal3,"MinimumVal3/D"); // Smallest voltage sample of 25um chamber signal
  outtree->Branch("DiffUncorr2",&DiffUncorr2,"DiffUncorr2/D"); // Difference of two consecutive Pedestal2 values
  outtree->Branch("DiffUncorr3",&DiffUncorr3,"DiffUncorr3/D"); // Difference of two consecutive Pedestal3 values
  outtree->Branch("Clusters2",&Clusters2); // Vector with index where cluster algorithm found a cluster, length = nClusters
  outtree->Branch("Clusters3",&Clusters3); // Vector with index where cluster algorithm found a cluster, length = nClusters

  outtree->Branch("StartIntegral2",&StartIntegral2,"StartIntegral2/I"); // Point from which the charge integration started.
  outtree->Branch("StartIntegral3",&StartIntegral3,"StartIntegral3/I"); // Point from which the charge integration started.
  outtree->Branch("FracCharge2",&FracCharge2,"FracCharge2/D"); // Fraction of total signal that was integrated.
  outtree->Branch("FracCharge3",&FracCharge3,"FracCharge3/D"); // Fraction of total signal that was integrated.
  outtree->Branch("RMS2",&prevRMS2,"RMS2/D"); // RMS of empty events in channel 2
  outtree->Branch("RMS3",&prevRMS3,"RMS3/D"); // RMS of empty events in channel 3
  outtree->Branch("Baseline2",&baseline2,"Baseline2/D"); // Average voltage of smoothed channel 2 signal for empty events.
  outtree->Branch("Baseline3",&baseline3,"Baseline3/D"); // Average voltage of smoothed channel 3 signal for empty events.

  outfile.cd();

  TGraph *gChan1[N_TRACEPLOTS];
  TGraph *gChan2[N_TRACEPLOTS];
  TGraph *gChan3[N_TRACEPLOTS];
  
  // Now we loop over all the entries
  //Int_t n_entries = C1->GetEntries();
  Int_t n_entries = 1000;
  for(Int_t i=0; i < n_entries; i++ )
    {
      // Re-initialize variables which are linked in the tree.
      SignalPresent2 = kFALSE;
      SignalPresent3 = kFALSE;
      rawTOF = -99;
      TOF = -99;
      nTOF = 0;
      Beta = 0;
      Charge2_pC = -9999;
      Charge3_pC = -9999;
      Pedestal1 = -9999;
      //Pedestal2_pC = -9999; These are only updated if nTOF == 0.
      //Pedestal3_pC = -9999;
      RawCharge2_pC = -9999;
      RawCharge3_pC = -9999;
      MinimumVal2 = 9999;
      MinimumVal3 = 9999;
      DiffUncorr2 = -9999;
      DiffUncorr3 = -9999;
      Clusters2.clear();
      Clusters3.clear();

      StartIntegral2 = -9999;
      StartIntegral3 = -9999;
      FracCharge2 = -1;
      FracCharge3 = -1;
      
      if(DEBUG == 1){cout << "Looping over trace plots...event " << i;}
      // Make trace plots.
      if( i < N_TRACEPLOTS )
	{
	  if(DEBUG == 1){cout << i;}
	  outfile.cd();
	  
	  C1->Draw("ampl_tv->fElements:Iteration$",TString::Format("Entry$==%i",i));
	  gChan1[i] = (TGraph*)gPad->GetPrimitive("Graph");
	  gChan1[i]->SetName(TString::Format("gChan1_%i",i));
	  gChan1[i]->SetTitle(TString::Format("Channel 1 Entry %i",i));
	  gChan1[i]->Write();

	  C2->Draw("ampl_tv->fElements:Iteration$",TString::Format("Entry$==%i",i));
	  gChan2[i] = (TGraph*)gPad->GetPrimitive("Graph");
	  gChan2[i]->SetName(TString::Format("gChan2_%i",i));
	  gChan2[i]->SetTitle(TString::Format("Channel 2 Entry %i",i));
	  gChan2[i]->Write();
	  
	  C3->Draw("ampl_tv->fElements:Iteration$",TString::Format("Entry$==%i",i));
	  gChan3[i] = (TGraph*)gPad->GetPrimitive("Graph");
	  gChan3[i]->SetName(TString::Format("gChan3_%i",i));
	  gChan3[i]->SetTitle(TString::Format("Channel 3 Entry %i",i));
	  gChan3[i]->Write();
	}
      if(DEBUG == 1){cout << ". Done. " << endl;}
      
      if(DEBUG == 1){cout << "Reading C1." << endl;}
      // Calculate the TOF from channel 1.
      f1.cd();
      C1->GetEntry(i);
      //..Get the pedestal from the first block of bins and modify 
      //  threshold accordingly
      C1_tv->GetSub(0,n_C1_pedestalframes-1,sub); 
      // GetSub(A,B) gets a subvector with elements A to B inclusive of the endpoints (hence the -1)
      Float_t pedestal = sub.Sum()/(1.0*n_C1_pedestalframes);
      Float_t mthresh = TOF_threshold+pedestal;
      Pedestal1 = pedestal;
      //..b0 to b3 are the four bins containing the negative transition across 
      //  the threshold
      Int_t b0 = 0, b1 = 0, b2 = 0, b3 = 0, db = 0;
      nTOF = 0;
      for(Int_t ib = 100; ib<20000; ib++ )
       	{
       	  if( (*C1_tv)[ib] < mthresh && (*C1_tv)[ib-1] > mthresh )
      	    {
      	      if( b0==0 ) { 
      	  	b0 = ib;
		nTOF++;
      	      } else if( b1==0 ) {
      	  	b1 = ib;
		nTOF++;
      	      } else if( b2==0 ) {
      	  	b2 = ib;
		nTOF++;
      	      } else if( b3==0 ) {
      	  	b3 = ib;
		nTOF++;
      	      }
	    }
	}
      db = b2+b3-b0-b1;
      rawTOF = db;
      TOF = (1.0*db/2.0*0.00000000005-0.00000006696)*1000000000.0; // TOF in nanosecond units with an offset of 66.96ns
      // Speed of particle (scintilator distance divided by time of flight) in units of c.
      Beta = (scintillator_distance/(1.0*db/2.0*0.00000000005-0.00000006696))/299800000.0; 
      
      if(DEBUG == 1){cout << "Reading C2." << endl;}
      f2.cd();
      C2->GetEntry(i);
      if(DEBUG == 1){cout << "Reading C3." <<  endl;}
      f3.cd();
      C3->GetEntry(i);

      smoothed2 = smoothed_tv(*C2_tv, n_frames_smooth);
      smoothed3 = smoothed_tv(*C3_tv, n_frames_smooth);

      newIntegralC2 = C2_tv->Sum();
      newIntegralC3 = C3_tv->Sum();
      
      // Here I apply a threshold algorithm on smoothed signals from C2 and C3.
      // I just take the first sample where the voltage deviates from the baseline
      // by more than 5 times the RMS (of the baseline).  

      for(Int_t j=0;!(SignalPresent2 && SignalPresent3) && j<N_SAMPLES;j++)
	{
	  if(!SignalPresent2 && ((smoothed2[j]-baseline2) < RMSMult*prevRMS2))
	    {
	      StartIntegral2 = max(j-prespace,0);// max is from std::algorithm, this makes the earliest integration time zero.
	      SignalPresent2 = kTRUE;
	    }
	  if(!SignalPresent3 && ((smoothed3[j]-baseline3) < RMSMult*prevRMS3))
	    {
	      StartIntegral3 = max(j-prespace,0);
	      SignalPresent3 = kTRUE;
	    }
	}
      if(StartIntegral2 == 0 || StartIntegral3 == 0)
	{
	  cout << "Skipping event " << i << ", the signal threshold algorithm triggered on the first sample." << endl;
	  continue;
	}
      
      if(nTOF == 0 && (SignalPresent2 || SignalPresent3))
       	{
      // 	  // If this is an uncorrelated event, but there is signal activity, skip the event.
       	  cout << "Event " << i << ", it has nTOF == 0 but signal activity.  Not skipping." << endl;
      // 	  continue;
       	}
      if(nTOF == 0)
	{
	  // This is evaluated if this an empty (or non-TOF) event.
	  // The baseline values are updated from the raw charge integral of this event.
	  // After this, the charge integrals are not done.

	  DiffUncorr2 = prevIntegralC2-newIntegralC2;
	  DiffUncorr3 = prevIntegralC3-newIntegralC3;

	  prevIntegralC2 = newIntegralC2;
	  prevIntegralC3 = newIntegralC3;
	  
	  baseline2 = smoothed2.Sum()/(1.0*N_SAMPLES);
	  baseline3 = smoothed3.Sum()/(1.0*N_SAMPLES);
	  // Here I calculate the RMS of the whole empty waveform for C2 and C3.
	  // The formula is the square root of the variance, and the variance is
	  // the mean of the squares minus the square of the mean.
	  prevRMS2 = TMath::Sqrt((smoothed2.Norm2Sqr()/(1.0*N_SAMPLES)) - TMath::Power(baseline2,2) );
	  prevRMS3 = TMath::Sqrt((smoothed3.Norm2Sqr()/(1.0*N_SAMPLES)) - TMath::Power(baseline3,2) );

	  Pedestal2_pC = -1.0*prevIntegralC2*sample_spacing/input_impedance;
	  Pedestal3_pC = -1.0*prevIntegralC3*sample_spacing/input_impedance;
	  if(DEBUG == 1){cout << "Re-measuring baselines at event " << i << endl;}
	}
      if(DEBUG == 1)
	{
	  cout << i << " channel 2 RMS: " << prevRMS2 << " baseline: " << baseline2 << " StartIntegral: " << StartIntegral2 << endl;
	  cout << i << " channel 3 RMS: " << prevRMS3 << " baseline: " << baseline3 << " StartIntegral: " << StartIntegral3 << endl;
	}
      if(prevIntegralC2 == 0 and prevIntegralC3 == 0)
	{
	  // This is evaluated if we haven't yet seen a single empty event.
	  // Since we don't have a usable baseline, we skip events until
	  // an empty one is seen.  This only wastes about 10 events.
	  cout << "Skipping event " << i << ", we don't have baselines measured yet." << endl;
	  continue;
	}
      if(nTOF == 4 && (!SignalPresent2 && !SignalPresent3))
	{
	  // This event has a good trigger, but has no signal activity in one or both chambers, so we skip it.
	  if(DEBUG == 1){cout << "Skipping event " << i << ", it has nTOF == 4, but no signal activity in one or both chambers." << endl;}
	  continue;
	}
      if(nTOF == 4 && SignalPresent2 && SignalPresent3)
	{
	  // This event has a good trigger and both chambers have signals.
	  // Here I do the charge integrations and cluster counting.
	  endintegral2 = min(N_SAMPLES,StartIntegral2+integral_duration); // The end of the integral is again not allowed to exceed the whole trace.
	  endintegral3 = min(N_SAMPLES,StartIntegral3+integral_duration);
	  // Some subvectors are created with the correct "duration".
	  C2_tv->GetSub(StartIntegral2,endintegral2-1,sub2);
	  C3_tv->GetSub(StartIntegral3,endintegral3-1,sub3);
	  // While prevIntegralC2/3 is a baseline to the subtracted from the rawcharge2/3 (see below), it is a baseline
	  // for an integral over the whole trace.  So here we scale the baseline by the fraction of samples which are
	  // being considered out of the whole trace.
	  sub_baseline2 = (endintegral2-StartIntegral2)/(1.0*N_SAMPLES)*prevIntegralC2;
	  sub_baseline3 = (endintegral3-StartIntegral3)/(1.0*N_SAMPLES)*prevIntegralC3;

	  // The actual "integrals"...
	  Charge2 = sub2.Sum();
	  Charge3 = sub3.Sum();

	  // Fraction of total charge integrated...
	  FracCharge2 = (Charge2-sub_baseline2)/(newIntegralC2-prevIntegralC2);
	  FracCharge3 = (Charge3-sub_baseline3)/(newIntegralC3-prevIntegralC3);

	  // Convert to picoCoulombs...
	  Charge2_pC = -1.0*(Charge2-sub_baseline2)*sample_spacing/input_impedance;
	  Charge3_pC = -1.0*(Charge3-sub_baseline3)*sample_spacing/input_impedance;

	  // "Raw" charge is the same as Charge but without baseline subtraction.
	  RawCharge2_pC = -1.0*Charge2*sample_spacing/input_impedance;
	  RawCharge3_pC = -1.0*Charge3*sample_spacing/input_impedance;

	  MinimumVal2 = C3_tv->Min();
	  MinimumVal3 = C3_tv->Min();

	  Clusters2 = clu_thresh_over_average(*C2_tv,cc_n_frames,cc_threshold);
	  Clusters3 = clu_thresh_over_average(*C3_tv,cc_n_frames,cc_threshold);
	}
      if(nTOF != 0 && nTOF != 4)
	{
	  cout << "Skipping event " << i << ", nTOF != (4,0)" << endl;
	  continue;
	}
 
      if(DEBUG == 1){cout << "Filling outtree." << endl;}
      outtree->Fill();
      
    } // Ends loop over n_entries
  
  if(DEBUG==1){cout << "Done looping over n_entries" << endl;}
  
  f1.Close();
  f2.Close();
  f3.Close();
  
  outfile.cd();
  gDirectory->Write();
  outfile.Close();
}

vector<Int_t> clu_thresh_over_average(TVectorD ampl_tv, Int_t n_frames, Float_t thresh)
{
  // An implementation of the algorithm described at
  //http://mailman.fe.infn.it/superbwiki/index.php/DCH_Cluster_Counting#Threshold_Over_Rolling_Average

  // ampl_tv is a TVectorD with the voltage samples for the DCH output.
  // n_frames is the number of frames to use in the rolling average.
  // thresh is a (generally positive) threshold used in the algorithm.

  // Form the N-frame rolling average.
  Int_t n_samples = ampl_tv.GetNoElements();
  // TVectorD avg(n_samples);
  // avg.Zero(); // Initialize a new TVectorD to contain the N-frame average.
  // Int_t n_averaging = n_frames;

  // for(Int_t i=0;i<n_samples;i++)
  //   {
  //     // The following line just takes care of the end of the sample, when there
  //     // are fewer samples remaining than n_frames.
  //     //n_averaging = (n_samples-i > n_frames) ? n_frames : n_samples-i;
  //     n_averaging = min(n_frames,n_samples-i);

  //     for(Int_t j=0;j<n_averaging;j++)
  // 	{
  // 	  avg[i] += ampl_tv[i+j];
  // 	}
  //     avg[i] /= n_averaging*1.0;
  //   }

  TVectorD avg = smoothed_tv(ampl_tv,n_frames);

  vector<Int_t> cluster_positions; // Vector to contain cluster positions.
  Int_t n_clusters = 0;
  for(Int_t i=2;i<n_samples;i++)
    {
      // If the non-averaged signal changed by more than thresh relative to the averaged
      // signal, then we count a cluster.
      if( (ampl_tv[i] - avg[i-1]) < thresh && (ampl_tv[i-1] - avg[i-2]) >= thresh )
	{
	  cluster_positions.push_back(i);
	  n_clusters++;
	}
    }
  return cluster_positions;
}

TVectorD smoothed_tv(TVectorD ampl_tv, Int_t n_frames)
{
  // Returns a smoothed version of the input vector.
  // ampl_tv is a TVectorD with the voltage samples for the DCH output.
  // n_frames is the number of frames to use in the rolling average.

  // Form the N-frame rolling average.
  Int_t n_samples = ampl_tv.GetNoElements();
  TVectorD smooth(n_samples);
  smooth.Zero(); // Initialize a new TVectorD to contain the N-frame average.

  for(Int_t i=0;i<=n_samples-n_frames;i++)
    {
      // The average is done "forwards", so we have to stop before i+n_frames = n_samples.

      // Add up the n_frames, then divide by n.
      for(Int_t j=0;j<n_frames;j++)
	{
	  smooth[i] += ampl_tv[i+j];
	}
      smooth[i] /= n_frames*1.0;
    }
  // For the remaining n_frames samples, we just use the last smoothed value.
  for(Int_t i=n_samples-n_frames+1;i<n_samples;i++)
    {
      smooth[i] = smooth[n_samples-n_frames];
    }

  return smooth;
}

Float_t get_scintilator_distance(TString path_to_dir)
{
  // This function tries to infer the run number from the path_to_dir which generally
  // ends with the run number, so the last 4 characters are taken as the run number.
  // Then different cases are tested to return the correct distance between the
  // scintillators used in that run.  If it fails, a warning is generated and the 
  // distance returned is the most common one of 3.884m.
  Int_t N = path_to_dir.Length();
  TString runstr = path_to_dir(N-4,N);
  Int_t runnum = runstr.Atoi();
  Float_t distance = 0;
  if(runnum <= 499 && runnum >= 300)
    {
      distance = 3.884;
    }
  else if (runnum == 500)
    {
      distance = 3.652;
    }
  else if (runnum == 501)
    {
      distance = 2.992;
    }
  else if (runnum >= 502 && runnum <= 516)
    {
      distance = 3.890;
    }
  else
    {
      cout << "Warning, did not find valid run number, using 3.884m scintilator spacing." << endl;
      distance = 3.884;
    }
  return distance;
}

int main(int argc, char * argv[])
{
  if(argc != 2)
    {
      cout << "Run the program with a path to a Run#### directory containing C1, C2, and C3 root files." << endl;
    }
  TString path_to_dir(argv[1]);
  standard_analysis(path_to_dir);
  return 0;
}