~jfcaron/+junk/TRIUMFBeamTest

// 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","path_to_outdir","iif",itime,cc_frames,cc_thresh)'
// Where path_to_dir is where are located the C1, C2 and C3.root files for the run that you wish to analyse.
// path_to_outdir is where you want the produced standard_analysis.root file.  "iif" is an example of a format
// string, with i and f representing integers and double floating point values in order.  The rest of the 
// variables must be of the right type as specified by the format string.
// 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 30 // Number of events for which to produce sample graphs.
#define DEBUG 0 // Set to 1 to see debug messages.
#define CC_ALGO 3
#define TEK 0 // Set to 0 for runs in the Aug-Sept 2012 test, set to 1 for December 2012.
#define CHANNEL3 0 // Set to 1 to require channel 3 to have signals along with channel 2

// This CPP section selects the right number of voltage samples.
#if TEK == 0
#define N_SAMPLES 20002 // The number of voltage samples in each oscilloscope trace. (LeCroy scope)
#define SAMPLE_WIDTH 50e-12 // 50 picoseconds per sample.
#define TOF_OFFSET 0.00000006696
#define SCALE 2
#elif TEK == 1
#define N_SAMPLES 4000 // Number of voltage samples per trace (Tektronix scope)
#define SCALE 4 // Ratio of this sample width to that for the LeCroy scope.
#define SAMPLE_WIDTH 200e-12 // 200ps per sample
#define TOF_OFFSET 6.6878e-8 // Determined by assuming that the electrons are moving at the speed of light.
#endif

// Which algorithm to use for cluster counting.  
//0: basic 2-parameter, 
//1: generalized 4-parameter,
//2: basic 2-parameter with timeout (3 parameter).
//3: second derivative
//4: generalized 4-parameter with timeout

#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"
#include "TTimeStamp.h"
#include "TSystem.h"
#include <cstdarg>

using namespace std;

// Include my own little library of algorithms.
#include "CC_algorithms.C"

// This is the main function which is run.
void standard_analysis(char * c_path_to_dir, char * c_path_to_out, char * c_argformat, ...)
{ 
  //gROOT->LoadMacro("CC_algorithms.C+");
  // Convert the c-string arguments to TStrings.
  TString path_to_dir(c_path_to_dir);
  TString path_to_out(c_path_to_out);
  if( path_to_out.IsNull() )
    {
      // If the second argument is an emtpy string "", default to putting
      // the output in the same directory as the input.
      path_to_out = path_to_dir;
    }
  TString argformat(c_argformat);

  // Get the runnum from the path.
  Int_t runnum = get_runnum(path_to_dir);

  // Keep a list of integer and double arguments.
  std::vector<Int_t> int_args;
  std::vector<Double_t> dbl_args;
  
  // Start building the outfilename.
  Int_t cc_algo = CC_ALGO;
  cout << "Using cluster-counting algorithm " << CC_ALGO << ", " << cc_algo << endl;
  TString outfilename = path_to_out+TString::Format("/standard_analyzed_Run%04d_CC%d",runnum,CC_ALGO);

  // This section uses macros from <cstdarg> to get the variable argument list.
  va_list args;
  va_start(args,c_argformat);
  Int_t int_arg = 0;
  Double_t dbl_arg = 0;
  for(Int_t i_arg = 0;i_arg < argformat.Length();i_arg++)
    {
      if( argformat[i_arg] == 'i')
	{
	  cout << i_arg << "th argument is integer: ";
	  int_arg = va_arg(args,Int_t);
	  cout << int_arg << endl;
	  int_args.push_back(int_arg);
	  outfilename.Append(TString::Format("_%d",int_arg));
	}
      else if (argformat[i_arg] == 'f')
	{
	  cout << i_arg << "th argument is double: ";
	  dbl_arg = va_arg(args,Double_t);
	  cout << dbl_arg << endl;
	  dbl_args.push_back(dbl_arg);
	  outfilename.Append(TString::Format("_%f",dbl_arg));
	}
    }
  outfilename.Append(TString(".root"));
  va_end(args);
  Int_t integration_time = int_args[0];

  // Get a unique filename based on the time
  //TTime theTime = gSystem->Now();
  TTimeStamp theTime = TTimeStamp();
  cout << "Processing started at " << theTime.AsString() << endl;
  cout << "Writing to outfile " << outfilename << endl;

  TFile outfile(outfilename,"RECREATE");

#if TEK == 0
  // Open files
  TFile f1(path_to_dir+"/C1.root");
  TFile f2(path_to_dir+"/C2.root");
  TFile f3(path_to_dir+"/C3.root");
  
  if (f1.IsZombie() || f2.IsZombie() || f3.IsZombie())
    {
      cout << "Error opening the CN.root files." << endl;
      exit(-1);
    }

 
  // 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;
    }
#elif TEK == 1
  // Open input file
  TFile f1(path_to_dir+"/"+TString::Format("Run%04d",runnum)+".root");
  if (f1.IsZombie())
    {
      cout << "Error opening the root file." << endl;
      exit(-1);
    }
  f1.cd();
  TTree* waves = (TTree*)gDirectory->Get("waves");
  Double_t C1_dbl[N_SAMPLES],C2_dbl[N_SAMPLES],C3_dbl[N_SAMPLES];
  waves->SetBranchAddress("ch1",&C1_dbl);
  waves->SetBranchAddress("ch2",&C2_dbl);
  waves->SetBranchAddress("ch3",&C3_dbl);
  TVectorD *C1_tv = 0;
  TVectorD *C2_tv = 0;
  TVectorD *C3_tv = 0;
#endif

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

  // Run-level information:
  // runnum was got higher before the files were opened.
  Bool_t channel3 = CHANNEL3;
  Int_t inverted_channels = get_inverted_channels(runnum);
  Int_t flipC2 = 1;
  Int_t flipC3 = 1;
  if (inverted_channels == 2)
    {
      flipC2 = -1;
    }
  else if(inverted_channels == 3)
    {
      flipC3 = -1;
    }

  // These are variables that are used during TOF calculations.
  Float_t scintillator_distance = get_scintilator_distance(runnum);
  const Int_t n_C1_pedestalframes = 1000/SCALE;
  Double_t Pedestal1 = 0; // Sum of first many samples of channel 1, updated at every event.
  Float_t TOF_threshold = -0.4; // For TOF calculation
  TVectorD sub(n_C1_pedestalframes); // Subvector for the TOF pedestal
  Int_t rawTOF = 0,nTOF = 0; // Unscaled TOF information, number of TOF threshold crossings
  Double_t TOF = 0,Beta = 0; // TOF information in nanoseconds, speed of particle in units of c.

  // These variables store the baselines and RMS of the empty events, used to set thresholds for the charge integrations.
  Float_t newIntegralC2 = 0, newIntegralC3 = 0; // This gets updated every event, it is the Sum() of all the voltage samples.
  Float_t prevIntegralC2 = 0; // This gets updated when we find an uncorrelated trigger event, it is the Sum() of all voltage samples.
  Float_t prevIntegralC3 = 0;
  Float_t chargeBaselineC2 = 0; // Holds the value for baseline subtraction of charge measurements.
  Float_t chargeBaselineC3 = 0;
  Double_t baseline2 = 0, baseline3 = 0; // In volts, the baseline of smoothed uncorrelated trigger signal
  Double_t prevRMS2 = 0, prevRMS3 = 0; // In volts, RMS of smoothed uncorrelated trigger signals
  Double_t Pedestal2_pC = 0,Pedestal3_pC = 0; // Integrated charge of uncorrelated trigger signals, in picoCoulombs.
  Double_t DiffUncorr2 = 0,DiffUncorr3 = 0; // Updated only at uncorrelated trigger events, difference between newIntegralC2 and prevIntegralC2.

  // These are constant values used in the signal threshold algorithm and charge integration
  Int_t prespace = 100/SCALE; // Frames to skip back from signal threshold for charge integral. 5 nanoseconds.
  const Int_t n_frames_smooth = 100/SCALE; // Number of frames of smoothing for signal threshold crossing, 5 nanoseconds.
  Int_t integral_duration = integration_time;
#if TEK == 0
  Float_t RMSMult = -5.0; // Fluctuation from baseline in units of RMS to trigger signal threshold.
  Float_t input_impedance = 75.0; // Ohms
#elif TEK == 1
  Float_t RMSMult = -3.0;
  Float_t input_impedance = 50.0; // Ohms
#endif
  Float_t sample_spacing = SAMPLE_WIDTH/1e-12; // In picoseconds.

  // These variables store values related to the signal threshold-crossing.
  Int_t StartIntegral2 = -99, StartIntegral3 = -99; // Where the signal threshold was crossed.
  Bool_t SignalPresent2 = kFALSE, SignalPresent3 = kFALSE; // Flags for threshold algorithm
  
  // These variables store values used for the charge integration step.
  TVectorD smoothed2(N_SAMPLES),smoothed3(N_SAMPLES); // These store the smoothed signals.
  Int_t endintegral2 = 0, endintegral3 = 0; // These store the endpoint of the integration
  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.
  Double_t Charge2_pC = 0,Charge3_pC = 0; // Scaled to pC
  Double_t RawCharge2_pC = 0,RawCharge3_pC = 0; // Charge in picoCoulombs without subtracting the baseline.
  Double_t FracCharge2 = 0, FracCharge3 = 0; // Fraction of total charge in the signal which was integrated into Charge2_pC.

  // These variables are used for the cluster-counting step.
  //Float_t cc_threshold = cc_thresh; // For cluster counting algorithm (in volts) Sam Dejong's "optimal" is 0.005
  //Int_t cc_n_frames = cc_frames; // For cluster counting algorithm (in frames) Sam Dejong's "optimal" is ~60
  vector<Int_t> Clusters2; // Array of indexes of clusters found by the algorithm.
  vector<Int_t> Clusters3;
  Double_t CluSep2 = 0; // Average separation of clusters found.
  Double_t CluSep3 = 0;

  // Additional calculated variables which will be linked to branches in outtree.
  Int_t Entry = 0; // Entry number, in case a weird event comes up, you can track it back to the original C# files or even .trc files.
  Double_t MinimumVal2 = 0,MinimumVal3 = 0; // Minimum voltage sample in a given trace, used to determine signal saturation in a run.

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

  // Run-level branches:
  outtree->Branch("runnum",&runnum,"runnum/I"); // Run number as an integer.
  outtree->Branch("AnalysisTime",&theTime,32000,0); // The time at which the standard_analysis function was called.
  outtree->Branch("IntegerArgs",&int_args); // Vector of integer arguments for cluster counting.
  outtree->Branch("DoubleArgs",&dbl_args); // Vector of double arguments for cluster counting.
  outtree->Branch("CC_ALGO",&cc_algo,"CC_ALGO/I"); // Code for cluster-counting algorithm.
  outtree->Branch("CHANNEL3",&channel3,"CHANNEL3/O"); // Whether channel 3 was considered in skipping events.

  // Entry-level branches:
  outtree->Branch("Entry",&Entry,"Entry/I"); // Entry number in the C1, C2, C3 trees.
  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 smoothed signal of empty events in channel 2
  outtree->Branch("RMS3",&prevRMS3,"RMS3/D"); // RMS of smoothed signal 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.
  outtree->Branch("CluSep2",&CluSep2,"CluSep2/D"); // Average separation of clusters found in channel 2.
  outtree->Branch("CluSep3",&CluSep3,"CluSep3/D"); // Average separation of clusters found in channel 3.

  outfile.cd();

  TGraph *gChan1[N_TRACEPLOTS];
  TGraph *gChan2[N_TRACEPLOTS];
  TGraph *gChan3[N_TRACEPLOTS];
  
  // Now we loop over all the entries
#if TEK == 0
  Int_t n_entries = C1->GetEntries();
#elif TEK == 1
  Int_t n_entries = waves->GetEntries();
#endif
  for(Int_t i=0; i < n_entries; i++ )
    {
      // Re-initialize some of the variables which are linked in the tree.
      Charge2_pC = -9999;
      Charge3_pC = -9999;
      RawCharge2_pC = -9999;
      RawCharge3_pC = -9999;
      MinimumVal2 = 9999;
      MinimumVal3 = 9999;
      DiffUncorr2 = -9999;
      DiffUncorr3 = -9999;
      FracCharge2 = -1;
      FracCharge3 = -1;

      Clusters2.clear();
      Clusters3.clear();

      CluSep2 = 0.0;
      CluSep3 = 0.0;

#if TEK == 0
      // Make trace plots.
      if( i < N_TRACEPLOTS )
      	{
	  if(DEBUG == 1){cout << "Looping over trace plots...event " << i << endl;}
      	  outfile.cd();

      	  C1->Draw("ampl_tv->fElements:Iteration$",TString::Format("Entry$==%i",i),"",1,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(TString::Format("ampl_tv->fElements*%i:Iteration$",flipC2),TString::Format("Entry$==%i",i),"",1,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(TString::Format("ampl_tv->fElements*%i:Iteration$",flipC3),TString::Format("Entry$==%i",i),"",1,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();
	}
#elif TEK == 1
      // Make trace plots.
      if( i < N_TRACEPLOTS )
      	{
	  if(DEBUG == 1){cout << "Looping over trace plots...event " << i << endl;}
      	  outfile.cd();

      	  waves->Draw("ch1->fElements:Iteration$",TString::Format("Entry$==%i",i),"",1,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();

      	  waves->Draw(TString::Format("ch2->fElements*%i:Iteration$",flipC2),TString::Format("Entry$==%i",i),"",1,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();
	  
      	  waves->Draw(TString::Format("ch3->fElements*%i:Iteration$",flipC3),TString::Format("Entry$==%i",i),"",1,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();
	}
#endif
#if TEK == 0
      if(DEBUG == 1){cout << "Reading C1." << endl;}
      // Calculate the TOF from channel 1.
      f1.cd();
      C1->GetEntry(i);
      if(DEBUG == 1){cout << "Reading C2." << endl;}
      f2.cd();
      C2->GetEntry(i);
      if(DEBUG == 1){cout << "Reading C3." <<  endl;}
      f3.cd();
      C3->GetEntry(i);
#elif TEK == 1
      f1.cd();
      waves->GetEntry(i);
      C1_tv = new TVectorD(N_SAMPLES,C1_dbl);
      C2_tv = new TVectorD(N_SAMPLES,C2_dbl);
      C3_tv = new TVectorD(N_SAMPLES,C3_dbl);
#endif

      if( inverted_channels == 2 )
	{
	  (*C2_tv)*= -1;
	}
      else if( inverted_channels == 3 )
	{
	  (*C3_tv)*= -1;
	}

      // We the pedestal from the first block of bins and set the actual 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/SCALE; ib<N_SAMPLES; 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*SAMPLE_WIDTH-TOF_OFFSET)*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*SAMPLE_WIDTH-TOF_OFFSET))/299800000.0; 

      if(nTOF != 0 && nTOF != 4)
	{
	  // These events have a mal-formed TOF value, so we skip them categorically.
	  // If they are very frequent, maybe something is wrong with the run?
	  cout << "Skipping event " << i << ", nTOF = " << nTOF << ", nTOF != (4,0)" << endl;
	  continue;
	}
      

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

      newIntegralC2 = C2_tv->Sum();
      newIntegralC3 = C3_tv->Sum();
      
      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.

	  if(DEBUG == 1){cout << "Re-measuring baselines at event " << i << endl;}
	  
	  DiffUncorr2 = prevIntegralC2-newIntegralC2;
	  DiffUncorr3 = prevIntegralC3-newIntegralC3;

	  // Here I set the baseline for charge integrals to the prevIntegral
	  // before it is re-updated in this empty event.  This allows the calculation
	  // of charge in empty events using the previous empty event as a baseline.
	  chargeBaselineC2 = prevIntegralC2;
	  chargeBaselineC3 = prevIntegralC3;

	  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;
	}
      else
	{
	  // If nTOF != 0, set charge baselines to the prevIntegralC2 value
	  chargeBaselineC2 = prevIntegralC2;
	  chargeBaselineC3 = prevIntegralC3;
	}
      if(DEBUG == 1)
	{
	  cout << i << " channel 2 RMS: " << prevRMS2 << 
	    " baseline: " << baseline2 << " StartIntegral: " << StartIntegral2 
	       << " chargeBaselineC2: " << chargeBaselineC2 << endl;
	  cout << i << " channel 3 RMS: " << prevRMS3 << 
	    " baseline: " << baseline3 << " StartIntegral: " << StartIntegral3 
	       << " chargeBaselineC2: " << chargeBaselineC2 << endl;
	}
      // This weird if construct is using the #defined CHANNEL3 constant which is
      // 0 or 1 depending on whether we wish CHANNEL3 to be important.  Here if 
      // CHANNEL3 is zero, the statement is just "and 1", if it is nonzero, then 
      // prevIntegralC3 is checked.  This is equivalent to the boolean statement
      // prevIntegralC2 == 0 and (prevIntegralC3 == 0 or !CHANNEL3)
      if(prevIntegralC2 == 0 and (CHANNEL3 ? (prevIntegralC3 == 0) : 1) )
	{
	  // 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.
	  // If more than ~10 events are skipped this way, there may be a
	  // problem with the run.
	  cout << "Skipping event " << i << ", we don't have baselines measured yet." << endl;
	  continue;
	}      

      // Here we apply a threshold algorithm on smoothed signals from C2 and C3.
      // We just take the first sample where the voltage deviates from the baseline
      // by more than RMSMult times the RMS (of the baseline).  RMSMult it set above.
      SignalPresent2 = kFALSE;
      SignalPresent3 = kFALSE;
      StartIntegral2 = -9999;
      StartIntegral3 = CHANNEL3 ? -9999 : 0;
      for(Int_t j=0;!(SignalPresent2 && SignalPresent3) && j<N_SAMPLES;j++) // The loop exits once both traces have a signal.
	{
	  if(!SignalPresent2 && ((smoothed2[j]-baseline2) < RMSMult*prevRMS2))
	    {
	      // max is from std::algorithm, this makes the earliest integration time zero.
	      StartIntegral2 = max(j-prespace,0);
	      SignalPresent2 = kTRUE;
	    }
	  if(!SignalPresent3 && ((smoothed3[j]-baseline3) < RMSMult*prevRMS3))
	    {
	      StartIntegral3 = max(j-prespace,0);
	      SignalPresent3 = kTRUE;
	    }
	}
      if(StartIntegral2 == 0 || (CHANNEL3 ? (StartIntegral3 == 0) : 0))
	{
	  // Sometime the baseline/RMS of the previous uncorrelated event is totally off
	  // from the correct baseline for the current event, so the threshold algorithm
	  // always triggers on the first sample that is tested (even if it's not sample #0).
	  // If this is the case, we have to throw away the event, since we don't have a good
	  // baseline.  If you see a lot of these, there may be a problem with the run.
	  cout << "Skipping event " << i << 
	    ", the signal threshold algorithm triggered on the first sample." << endl;
	  continue;
	}
      
      if(nTOF == 0 && (SignalPresent2 || (CHANNEL3 ? SignalPresent3 : 0)))
       	{
       	  // Sometimes we have an uncorrelated trigger event that nonetheless has signal activity
	  // (according to the threshold algorithm).  We could throw away these events, but actually
	  // most of these tend to be valid uncorrelated trigger events with an accidental signal
	  // threshold crossing, rather than real signals contaminating our empty events.  If we skip
	  // these events, it tends to throw off the baseline calculation.  In your own code you may
	  // wish to skip these events after all, but it's not immediately obvious that it helps.
	  // As before, if there are a lot of these events, maybe there is a bigger problem with the run.
       	  cout << "Event " << i << ", it has nTOF == 0 but signal activity.  Not skipping. Channels 2,3: " 
	       << SignalPresent2 << SignalPresent3 << endl;
       	}

      if(nTOF == 4 && (!SignalPresent2 or (CHANNEL3 ? !SignalPresent3 : 0)))
	{
	  // This event has a good trigger, but has no signal activity in one or both chambers, so we skip it.
	  // This is actually relatively common, so the cout is not performed normally.  It happens when 
	  // particles trigger the TOF system but either miss the chamber or leave no discernible signal in the
	  // chambers.  
	  if(DEBUG == 1){cout << "Skipping event " << i << ", it has nTOF == 4, but no signal activity in one or both chambers." << endl;}
	  //cout << "Skipping event " << i << ", it has nTOF == 4, but no signal activity in one or both chambers." << endl;
	  continue;
	}
      if(nTOF == 0 && !(SignalPresent2 || SignalPresent3))
	{
	  // This is an empty event that did not trigger the signal threshold algorithm.
	  StartIntegral2 = 2000/SCALE;
	  StartIntegral3 = 2000/SCALE;
	}
      if(StartIntegral2 >= 0 && (CHANNEL3 ? (StartIntegral3 >= 0) : 1))
	{
	  // 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".
	  TVectorD sub2,sub3;
	  C2_tv->GetSub(StartIntegral2,endintegral2-1,sub2);
	  C3_tv->GetSub(StartIntegral3,endintegral3-1,sub3);

	  // While chargeBaselineC2/3 is a baseline to be 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)*chargeBaselineC2;
	  sub_baseline3 = (endintegral3-StartIntegral3)/(1.0*N_SAMPLES)*chargeBaselineC3;

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

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

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

	  // Minimum value of each signal, for saturation study.
	  MinimumVal2 = C3_tv->Min();
	  MinimumVal3 = C3_tv->Min();

	  // The cluster-counting algorithm is implemented as a separate function.
	  // To add new cluster counting algorithms, this part needs to be edited.
	  switch ( CC_ALGO )
	    {
	    case 0:
	      // Basic 2-parameter algorithm
	      // Argument 1: the digitized signal.
	      // Argument 2: number of frames over which to smooth: int_args[1]
	      // Argument 3: threshold to apply in mV: dlb_args[0]
	      Clusters2 = clu_thresh_over_average(*C2_tv,int_args[1],dbl_args[0]);
	      Clusters3 = clu_thresh_over_average(*C3_tv,int_args[1],dbl_args[0]);
	      CluSep2 = Cellwise_Separation(Clusters2);
	      CluSep3 = Cellwise_Separation(Clusters3);
	      break;
	    case 1:
	      // More general 4-parameter algorithm
	      // Argument 1: the digitized signal.
	      // Argument 2: number of frames for smoothing "signal": int_args[1]
	      // Argument 3: number of frames for smoothing "baseline": int_args[2]
	      // Argument 4: number of frames to delay the difference between signal and baseline: int_args[3]
	      // Argument 5: threshold to apply to difference in mV: dbl_args[0]
	      Clusters2 = clu_hardware_derivative(*C2_tv,int_args[1],int_args[2],int_args[3],dbl_args[0]);
	      Clusters3 = clu_hardware_derivative(*C3_tv,int_args[1],int_args[2],int_args[3],dbl_args[0]);
	      CluSep2 = Cellwise_Separation(Clusters2);
	      CluSep3 = Cellwise_Separation(Clusters3);
	      break;
	    case 2:
	      // Same as CC_ALGO == 0, but with a timeout for removing fake clusters.
	      // Argument 1: the digitized signal.
	      // Argument 2: number of frames over which to smooth: int_args[1]
	      // Argument 3: number of frames that clusters need to stay below their initial voltage.
	      // Argument 4: threshold to apply in mV: dlb_args[0]
	      Clusters2 = timeout_booster(*C2_tv, clu_thresh_over_average(*C2_tv,int_args[1],dbl_args[0]), int_args[2]);
	      Clusters3 = timeout_booster(*C3_tv, clu_thresh_over_average(*C3_tv,int_args[1],dbl_args[0]), int_args[2]);
	      CluSep2 = Cellwise_Separation(Clusters2);
	      CluSep3 = Cellwise_Separation(Clusters3);
	      break;
	    case 3:
	      // This is the second-derivative algorithm.  It has two parameters.
	      // Argument 1: digitized signal
	      // Argument 2: number of frames over which to smooth, using the smoothed_reduce algorithm.
	      // Argument 3: threshold to apply to the 2nd derivative, in V/(50ps)**2.
	      Clusters2 = clu_second_derivative(*C2_tv,int_args[1],dbl_args[0]);
	      Clusters3 = clu_second_derivative(*C3_tv,int_args[1],dbl_args[0]);
	      CluSep2 = Cellwise_Separation(Clusters2);
	      CluSep3 = Cellwise_Separation(Clusters3);
	      break;
	    case 4:
	      // This is the same as CC_ALGO == 1, but with the timeout booster applied.
	      // int_args[0] : unused, it is for dE/dx integration
	      // int_args[1] : number of frames for smoothing "signal"
	      // int_args[2] : number of frames for smoothing "baseline"
	      // int_args[3] : number of frames to delay the difference between signal and baseline
	      // int_args[4] : timeout to use for booster
	      // dbl_args[0] : threshold to apply.
	      Clusters2 = timeout_booster(*C2_tv, clu_hardware_derivative(*C2_tv,int_args[1],int_args[2],int_args[3],dbl_args[0]),
					  int_args[4]);
	      Clusters3 = timeout_booster(*C3_tv, clu_hardware_derivative(*C3_tv,int_args[1],int_args[2],int_args[3],dbl_args[0]),
					  int_args[4]);
	      CluSep2 = Cellwise_Separation(Clusters2);
	      CluSep3 = Cellwise_Separation(Clusters3);
	      break;
	    case 5:
	      // This is the two-pass second derivative method.  It has the same parameters
	      // as CC_ALGO 3:
	      // Argument 1: digitized signal
	      // Argument 2: number of frames over which to smooth, using the smoothed_reduce algorithm.
	      // Argument 3: threshold to apply to the 2nd derivative, in V/(50ps)**2.
	      Clusters2 = clu_second_derivative(*C2_tv,int_args[1],dbl_args[0]);
	      Clusters3 = clu_second_derivative(*C3_tv,int_args[1],dbl_args[0]);
	      
	      vector<Int_t> BonusClusters2 = clu_second_derivative(*C2_tv,int_args[1],dbl_args[0],int_args[1]/2);
	      vector<Int_t> BonusClusters3 = clu_second_derivative(*C3_tv,int_args[1],dbl_args[0],int_args[1]/2);

	      CluSep2 = Combine_Separations(Cellwise_Separation(Clusters2), Clusters2.size(), 
					    Cellwise_Separation(BonusClusters2), BonusClusters2.size());
	      CluSep3 = Combine_Separations(Cellwise_Separation(Clusters3), Clusters3.size(), 
					    Cellwise_Separation(BonusClusters3), BonusClusters3.size());

	      Clusters2.insert(Clusters2.begin(),BonusClusters2.begin(),BonusClusters2.end());
	      Clusters3.insert(Clusters3.begin(),BonusClusters3.begin(),BonusClusters3.end());

	      sort(Clusters2.begin(), Clusters2.end());
	      sort(Clusters3.begin(), Clusters3.end());
	      break;
	    }
	  
	}
      
      Entry = i;

      if(DEBUG == 1){cout << "Filling outtree." << endl;}
      outtree->Fill();
#if TEK == 1
      delete C1_tv;
      delete C2_tv;
      delete C3_tv;
#endif
      
    } // Ends loop over n_entries
  
  if(DEBUG==1){cout << "Done looping over n_entries" << endl;}
  
  f1.Close();
#if TEK == 0
  f2.Close();
  f3.Close();
#endif  
  
  outfile.cd();
  gDirectory->Write();
  outfile.Close();

  cout << "Finished writing to file " << outfilename << endl;
}


// This main function exists so that the code can be compiled with gcc.  ROOT documentation warns
// to NOT define a main function if compiling with ACLiC, so it is commented out.

//The following compilation command should work (from bash):
// g++ -Wall -O3 -o standard_analysis standard_analysis.C -I`root-config --incdir` `root-config --glibs`
// The optimization level seems to be irrelevant beyond -O2, but -O0 is very slow and -O1 also slow.
// Note that it is still faster to run through ACLiC/ROOT rather than compile and run separately, but
// it may be easier to run it on a distributed system this way.
// int main(int argc, char * argv[])
// {
//   TString path_to_dir(argv[1]);
//   TString argformat("iid");
//   standard_analysis(path_to_dir, argformat, integration_time, cc_frames, cc_thresh);
//   return 0;
// }

// This is an overloaded function with the same name as the primary function, but has no arguments or content.
// It allows one to compile the function with root -q -b -l standard_analysis.C++ without arguments.
//void standard_analysis() {}