From e0d9b2a9b7c3b06b4f3a536fc6aab38a54f4f9d5 Mon Sep 17 00:00:00 2001
From: Eric Pooser <pooser@jlab.org>
Date: Mon, 20 Mar 2017 15:55:01 -0400
Subject: [PATCH] Complete rewrite of Aerogel and Cerenkov detector classes

New parameters are added.  hallc_replay of 3/31/2017 or later is needed.

As of now, the size of the TClonesArray objects are hard-coded in the
respective header files. This is not okay in the long run. Ideally,
the size of the TClonesArrays should be defined by the parameter which
dictates how many readout detectors exist for the respective
detector. I tried doing this in the Init method however, podd
complained regarding the size of the TClonesArrays in the DEF-file
histograms. I was unsuccessful in my first attempt. This should be
addressed in the future even though things work for now.
---
 src/THcAerogel.cxx   | 1081 +++++++++++++++++++++++++-----------------
 src/THcAerogel.h     |  197 +++++---
 src/THcCherenkov.cxx |  546 +++++++++++----------
 src/THcCherenkov.h   |  134 +++---
 4 files changed, 1124 insertions(+), 834 deletions(-)

diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx
index 11c74d9..2b0cffc 100644
--- a/src/THcAerogel.cxx
+++ b/src/THcAerogel.cxx
@@ -1,11 +1,10 @@
 /** \class THcAerogel
     \ingroup Detectors
 
-Class for an Aerogel detector consisting of pairs of PMT's
-attached to a diffuser box
-Will have a fixed number of pairs, but need to later make this
-configurable.
-
+    Class for an Aerogel detector consisting of pairs of PMT's
+    attached to a diffuser box
+    Will have a fixed number of pairs, but need to later make this
+    configurable.
 */
 
 #include "THcAerogel.h"
@@ -24,7 +23,6 @@ configurable.
 #include "THaTrack.h"
 #include "TClonesArray.h"
 #include "TMath.h"
-
 #include "THaTrackProj.h"
 #include "THcRawAdcHit.h"
 
@@ -40,33 +38,55 @@ THcAerogel::THcAerogel( const char* name, const char* description,
                         THaApparatus* apparatus ) :
   THaNonTrackingDetector(name,description,apparatus)
 {
-  // Normal constructor with name and description
-  fPosTDCHits = new TClonesArray("THcSignalHit",16);
-  fNegTDCHits = new TClonesArray("THcSignalHit",16);
-  fPosADCHits = new TClonesArray("THcSignalHit",16);
-  fNegADCHits = new TClonesArray("THcSignalHit",16);
-
-  frPosAdcPedRaw = new TClonesArray("THcSignalHit", 16);
-  frPosAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
-  frPosAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
-  frPosAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
-
-  frPosAdcPed = new TClonesArray("THcSignalHit", 16);
-  frPosAdcPulseInt = new TClonesArray("THcSignalHit", 16);
-  frPosAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
 
-  frNegAdcPedRaw = new TClonesArray("THcSignalHit", 16);
-  frNegAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
-  frNegAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
-  frNegAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
-
-  frNegAdcPed = new TClonesArray("THcSignalHit", 16);
-  frNegAdcPulseInt = new TClonesArray("THcSignalHit", 16);
-  frNegAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
+  // Normal constructor with name and description
+  frPosAdcPedRaw       = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPulseIntRaw  = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPulseAmpRaw  = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPulseTimeRaw = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPed          = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPulseInt     = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frPosAdcPulseAmp     = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  frNegAdcPedRaw       = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPulseIntRaw  = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPulseAmpRaw  = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPulseTimeRaw = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPed          = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPulseInt     = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  frNegAdcPulseAmp     = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+  fPosAdcErrorFlag     = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  fNegAdcErrorFlag     = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+
+  fNumPosAdcHits         = vector<Int_t>    (MaxNumPosAeroPmt, 0.0);
+  fNumGoodPosAdcHits     = vector<Int_t>    (MaxNumPosAeroPmt, 0.0);
+  fNumNegAdcHits         = vector<Int_t>    (MaxNumNegAeroPmt, 0.0);
+  fNumGoodNegAdcHits     = vector<Int_t>    (MaxNumNegAeroPmt, 0.0);
+  fNumTracksMatched      = vector<Int_t>    (MaxNumPosAeroPmt, 0.0);
+  fNumTracksFired        = vector<Int_t>    (MaxNumPosAeroPmt, 0.0);
+  fPosNpe                = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fNegNpe                = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodPosAdcPed         = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodPosAdcPulseInt    = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodPosAdcPulseIntRaw = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodPosAdcPulseAmp    = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodPosAdcPulseTime   = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fGoodNegAdcPed         = vector<Double_t> (MaxNumNegAeroPmt, 0.0);
+  fGoodNegAdcPulseInt    = vector<Double_t> (MaxNumNegAeroPmt, 0.0);
+  fGoodNegAdcPulseIntRaw = vector<Double_t> (MaxNumNegAeroPmt, 0.0);
+  fGoodNegAdcPulseAmp    = vector<Double_t> (MaxNumNegAeroPmt, 0.0);
+  fGoodNegAdcPulseTime   = vector<Double_t> (MaxNumNegAeroPmt, 0.0);
+
+  // 6 GeV variables
+  fPosTDCHits = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*16);
+  fNegTDCHits = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*16);
+  fPosADCHits = new TClonesArray("THcSignalHit", MaxNumPosAeroPmt*MaxNumAdcPulse);
+  fNegADCHits = new TClonesArray("THcSignalHit", MaxNumNegAeroPmt*MaxNumAdcPulse);
+
+  fPosNpeSixGev = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
+  fNegNpeSixGev = vector<Double_t> (MaxNumPosAeroPmt, 0.0);
 
   InitArrays();
 
-//  fTrackProj = new TClonesArray( "THaTrackProj", 5 );
 }
 
 //_____________________________________________________________________________
@@ -74,29 +94,29 @@ THcAerogel::THcAerogel( ) :
   THaNonTrackingDetector()
 {
   // Constructor
+  frPosAdcPedRaw       = NULL;
+  frPosAdcPulseIntRaw  = NULL;
+  frPosAdcPulseAmpRaw  = NULL;
+  frPosAdcPulseTimeRaw = NULL;
+  frPosAdcPed          = NULL;
+  frPosAdcPulseInt     = NULL;
+  frPosAdcPulseAmp     = NULL;
+  frNegAdcPedRaw       = NULL;
+  frNegAdcPulseIntRaw  = NULL;
+  frNegAdcPulseAmpRaw  = NULL;
+  frNegAdcPulseTimeRaw = NULL;
+  frNegAdcPed          = NULL;
+  frNegAdcPulseInt     = NULL;
+  frNegAdcPulseAmp     = NULL;
+  fPosAdcErrorFlag     = NULL;
+  fNegAdcErrorFlag     = NULL;
+
+  // 6 GeV variables
   fPosTDCHits = NULL;
   fNegTDCHits = NULL;
   fPosADCHits = NULL;
   fNegADCHits = NULL;
 
-  frPosAdcPedRaw = NULL;
-  frPosAdcPulseIntRaw = NULL;
-  frPosAdcPulseAmpRaw = NULL;
-  frPosAdcPulseTimeRaw = NULL;
-
-  frPosAdcPed = NULL;
-  frPosAdcPulseInt = NULL;
-  frPosAdcPulseAmp = NULL;
-
-  frNegAdcPedRaw = NULL;
-  frNegAdcPulseIntRaw = NULL;
-  frNegAdcPulseAmpRaw = NULL;
-  frNegAdcPulseTimeRaw = NULL;
-
-  frNegAdcPed = NULL;
-  frNegAdcPulseInt = NULL;
-  frNegAdcPulseAmp = NULL;
-
   InitArrays();
 
 }
@@ -105,98 +125,99 @@ THcAerogel::THcAerogel( ) :
 THcAerogel::~THcAerogel()
 {
   // Destructor
-  DeleteArrays();
+  delete frPosAdcPedRaw;       frPosAdcPedRaw       = NULL;
+  delete frPosAdcPulseIntRaw;  frPosAdcPulseIntRaw  = NULL;
+  delete frPosAdcPulseAmpRaw;  frPosAdcPulseAmpRaw  = NULL;
+  delete frPosAdcPulseTimeRaw; frPosAdcPulseTimeRaw = NULL;
+  delete frPosAdcPed;          frPosAdcPed          = NULL;
+  delete frPosAdcPulseInt;     frPosAdcPulseInt     = NULL;
+  delete frPosAdcPulseAmp;     frPosAdcPulseAmp     = NULL;
+  delete frNegAdcPedRaw;       frNegAdcPedRaw       = NULL;
+  delete frNegAdcPulseIntRaw;  frNegAdcPulseIntRaw  = NULL;
+  delete frNegAdcPulseAmpRaw;  frNegAdcPulseAmpRaw  = NULL;
+  delete frNegAdcPulseTimeRaw; frNegAdcPulseTimeRaw = NULL;
+  delete frNegAdcPed;          frNegAdcPed          = NULL;
+  delete frNegAdcPulseInt;     frNegAdcPulseInt     = NULL;
+  delete frNegAdcPulseAmp;     frNegAdcPulseAmp     = NULL;
+  delete fPosAdcErrorFlag;     fPosAdcErrorFlag     = NULL;
+  delete fNegAdcErrorFlag;     fNegAdcErrorFlag     = NULL;
 
+  // 6 GeV variables
   delete fPosTDCHits; fPosTDCHits = NULL;
   delete fNegTDCHits; fNegTDCHits = NULL;
   delete fPosADCHits; fPosADCHits = NULL;
   delete fNegADCHits; fNegADCHits = NULL;
 
-  delete frPosAdcPedRaw; frPosAdcPedRaw = NULL;
-  delete frPosAdcPulseIntRaw; frPosAdcPulseIntRaw = NULL;
-  delete frPosAdcPulseAmpRaw; frPosAdcPulseAmpRaw = NULL;
-  delete frPosAdcPulseTimeRaw; frPosAdcPulseTimeRaw = NULL;
-
-  delete frPosAdcPed; frPosAdcPed = NULL;
-  delete frPosAdcPulseInt; frPosAdcPulseInt = NULL;
-  delete frPosAdcPulseAmp; frPosAdcPulseAmp = NULL;
-
-  delete frNegAdcPedRaw; frNegAdcPedRaw = NULL;
-  delete frNegAdcPulseIntRaw; frNegAdcPulseIntRaw = NULL;
-  delete frNegAdcPulseAmpRaw; frNegAdcPulseAmpRaw = NULL;
-  delete frNegAdcPulseTimeRaw; frNegAdcPulseTimeRaw = NULL;
+  DeleteArrays();
 
-  delete frNegAdcPed; frNegAdcPed = NULL;
-  delete frNegAdcPulseInt; frNegAdcPulseInt = NULL;
-  delete frNegAdcPulseAmp; frNegAdcPulseAmp = NULL;
 }
 
 //_____________________________________________________________________________
 void THcAerogel::InitArrays()
 {
-  fA_Pos = NULL;
-  fA_Neg = NULL;
-  fA_Pos_p = NULL;
-  fA_Neg_p = NULL;
-  fT_Pos = NULL;
-  fT_Neg = NULL;
   fPosGain = NULL;
   fNegGain = NULL;
+
+  // 6 GeV variables
+  fA_Pos       = NULL;
+  fA_Neg       = NULL;
+  fA_Pos_p     = NULL;
+  fA_Neg_p     = NULL;
+  fT_Pos       = NULL;
+  fT_Neg       = NULL;
   fPosPedLimit = NULL;
   fNegPedLimit = NULL;
-  fPosPedMean = NULL;
-  fNegPedMean = NULL;
-  fPosNpe = NULL;
-  fNegNpe = NULL;
-  fPosPedSum = NULL;
-  fPosPedSum2 = NULL;
+  fPosPedMean  = NULL;
+  fNegPedMean  = NULL;
+  fPosPedSum   = NULL;
+  fPosPedSum2  = NULL;
   fPosPedCount = NULL;
-  fNegPedSum = NULL;
-  fNegPedSum2 = NULL;
+  fNegPedSum   = NULL;
+  fNegPedSum2  = NULL;
   fNegPedCount = NULL;
-  fPosPed = NULL;
-  fPosSig = NULL;
-  fPosThresh = NULL;
-  fNegPed = NULL;
-  fNegSig = NULL;
-  fNegThresh = NULL;
+  fPosPed      = NULL;
+  fPosSig      = NULL;
+  fPosThresh   = NULL;
+  fNegPed      = NULL;
+  fNegSig      = NULL;
+  fNegThresh   = NULL;
 }
 //_____________________________________________________________________________
 void THcAerogel::DeleteArrays()
 {
-  delete [] fA_Pos; fA_Pos = NULL;
-  delete [] fA_Neg; fA_Neg = NULL;
-  delete [] fA_Pos_p; fA_Pos_p = NULL;
-  delete [] fA_Neg_p; fA_Neg_p = NULL;
-  delete [] fT_Pos; fT_Pos = NULL;
-  delete [] fT_Neg; fT_Neg = NULL;
   delete [] fPosGain; fPosGain = NULL;
   delete [] fNegGain; fNegGain = NULL;
+
+  // 6 GeV variables
+  delete [] fA_Pos;       fA_Pos       = NULL;
+  delete [] fA_Neg;       fA_Neg       = NULL;
+  delete [] fA_Pos_p;     fA_Pos_p     = NULL;
+  delete [] fA_Neg_p;     fA_Neg_p     = NULL;
+  delete [] fT_Pos;       fT_Pos       = NULL;
+  delete [] fT_Neg;       fT_Neg       = NULL;
   delete [] fPosPedLimit; fPosPedLimit = NULL;
   delete [] fNegPedLimit; fNegPedLimit = NULL;
-  delete [] fPosPedMean; fPosPedMean = NULL;
-  delete [] fNegPedMean; fNegPedMean = NULL;
-  delete [] fPosNpe; fPosNpe = NULL;
-  delete [] fNegNpe; fNegNpe = NULL;
-  delete [] fPosPedSum; fPosPedSum = NULL;
-  delete [] fPosPedSum2; fPosPedSum2 = NULL;
+  delete [] fPosPedMean;  fPosPedMean  = NULL;
+  delete [] fNegPedMean;  fNegPedMean  = NULL;
+  delete [] fPosPedSum;   fPosPedSum   = NULL;
+  delete [] fPosPedSum2;  fPosPedSum2  = NULL;
   delete [] fPosPedCount; fPosPedCount = NULL;
-  delete [] fNegPedSum; fNegPedSum = NULL;
-  delete [] fNegPedSum2; fNegPedSum2 = NULL;
+  delete [] fNegPedSum;   fNegPedSum   = NULL;
+  delete [] fNegPedSum2;  fNegPedSum2  = NULL;
   delete [] fNegPedCount; fNegPedCount = NULL;
-  delete [] fPosPed; fPosPed = NULL;
-  delete [] fPosSig; fPosSig = NULL;
-  delete [] fPosThresh; fPosThresh = NULL;
-  delete [] fNegPed; fNegPed = NULL;
-  delete [] fNegSig; fNegSig = NULL;
-  delete [] fNegThresh; fNegThresh = NULL;
+  delete [] fPosPed;      fPosPed      = NULL;
+  delete [] fPosSig;      fPosSig      = NULL;
+  delete [] fPosThresh;   fPosThresh   = NULL;
+  delete [] fNegPed;      fNegPed      = NULL;
+  delete [] fNegSig;      fNegSig      = NULL;
+  delete [] fNegThresh;   fNegThresh   = NULL;
 }
 
 //_____________________________________________________________________________
 THaAnalysisObject::EStatus THcAerogel::Init( const TDatime& date )
 {
 
-  cout << "THcAerogel::Init " << GetName() << endl;
+  cout << "THcAerogel::Init for: " << GetName() << endl;
 
   char EngineDID[] = "xAERO";
   EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
@@ -223,55 +244,94 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
   // This function is called by THaDetectorBase::Init() once at the beginning
   // of the analysis.
 
-  cout << "THcAerogel::ReadDatabase " << GetName() << endl;
+  cout << "THcAerogel::ReadDatabase for: " << GetName() << endl;
 
   char prefix[2];
-
   prefix[0]=tolower(GetApparatus()->GetName()[0]);
   prefix[1]='\0';
 
-  fNelem = 8;			// Default if not defined
-  Bool_t optional=true ;
+  fNRegions = 1;   // Default if not set in parameter file
+
   DBRequest listextra[]={
-    {"aero_num_pairs", &fNelem, kInt,0,optional},
+    {"aero_num_pairs", &fNelem, kInt},
     {0}
   };
-  gHcParms->LoadParmValues((DBRequest*)&listextra,prefix);
 
-  fA_Pos = new Float_t[fNelem];
-  fA_Neg = new Float_t[fNelem];
-  fA_Pos_p = new Float_t[fNelem];
-  fA_Neg_p = new Float_t[fNelem];
-  fT_Pos = new Float_t[fNelem];
-  fT_Neg = new Float_t[fNelem];
+  gHcParms->LoadParmValues((DBRequest*)&listextra, prefix);
+
+  Bool_t optional = true ;
+
+  cout << "Number of " << GetApparatus()->GetName() << "."
+       << GetName() << " PMT pairs defined = " << fNelem << endl;
 
   fPosGain = new Double_t[fNelem];
   fNegGain = new Double_t[fNelem];
+
+  // 6 GeV variables
+  fTdcOffset   = 0; // Offset to make reference time subtracted times positve
   fPosPedLimit = new Int_t[fNelem];
   fNegPedLimit = new Int_t[fNelem];
-  fPosPedMean = new Double_t[fNelem];
-  fNegPedMean = new Double_t[fNelem];
-
-  fTdcOffset = 0;		// Offset to make reference time subtracted times positve
+  fA_Pos       = new Float_t[fNelem];
+  fA_Neg       = new Float_t[fNelem];
+  fA_Pos_p     = new Float_t[fNelem];
+  fA_Neg_p     = new Float_t[fNelem];
+  fT_Pos       = new Float_t[fNelem];
+  fT_Neg       = new Float_t[fNelem];
+  fPosPedMean  = new Double_t[fNelem];
+  fNegPedMean  = new Double_t[fNelem];
 
   // Create arrays to hold pedestal results
-  InitializePedestals();
+  if (fSixGevData) InitializePedestals();
+
+  // Region parameters
+  fRegionsValueMax = fNRegions * 8;
+  fRegionValue     = new Double_t[fRegionsValueMax];
 
   DBRequest list[]={
-    {"aero_pos_gain", fPosGain, kDouble, (UInt_t) fNelem},
-    {"aero_neg_gain", fNegGain, kDouble, (UInt_t) fNelem},
-    {"aero_pos_ped_limit", fPosPedLimit, kInt, (UInt_t) fNelem},
-    {"aero_neg_ped_limit", fNegPedLimit, kInt, (UInt_t) fNelem},
-    {"aero_pos_ped_mean", fPosPedMean, kDouble, (UInt_t) fNelem,optional},
-    {"aero_neg_ped_mean", fNegPedMean, kDouble, (UInt_t) fNelem,optional},
-    {"aero_tdc_offset", &fTdcOffset, kInt, 0, optional},
-    {"aero_min_peds", &fMinPeds, kInt, 0, optional},
+    {"aero_num_regions",      &fNRegions,         kInt},
+    {"aero_red_chi2_min",     &fRedChi2Min,       kDouble},
+    {"aero_red_chi2_max",     &fRedChi2Max,       kDouble},
+    {"aero_beta_min",         &fBetaMin,          kDouble},
+    {"aero_beta_max",         &fBetaMax,          kDouble},
+    {"aero_enorm_min",        &fENormMin,         kDouble},
+    {"aero_enorm_max",        &fENormMax,         kDouble},
+    {"aero_diff_box_zpos",    &fDiffBoxZPos,      kDouble},
+    {"aero_npe_thresh",       &fNpeThresh,        kDouble},
+    {"aero_adcTimeWindowMin", &fAdcTimeWindowMin, kDouble},
+    {"aero_adcTimeWindowMax", &fAdcTimeWindowMax, kDouble},
+    {"aero_debug_adc",        &fDebugAdc,         kInt,    0, 1},
+    {"aero_six_gev_data",     &fSixGevData,       kInt,    0, 1},
+    {"aero_pos_gain",         fPosGain,           kDouble, (UInt_t) fNelem},
+    {"aero_neg_gain",         fNegGain,           kDouble, (UInt_t) fNelem},
+    {"aero_pos_ped_limit",    fPosPedLimit,       kInt,    (UInt_t) fNelem, optional},
+    {"aero_neg_ped_limit",    fNegPedLimit,       kInt,    (UInt_t) fNelem, optional},
+    {"aero_pos_ped_mean",     fPosPedMean,        kDouble, (UInt_t) fNelem, optional},
+    {"aero_neg_ped_mean",     fNegPedMean,        kDouble, (UInt_t) fNelem, optional},
+    {"aero_tdc_offset",       &fTdcOffset,        kInt,    0,               optional},
+    {"aero_min_peds",         &fMinPeds,          kInt,    0,               optional},
+    {"aero_region",           &fRegionValue[0],   kDouble, (UInt_t) fRegionsValueMax},
+
     {0}
   };
-  gHcParms->LoadParmValues((DBRequest*)&list,prefix);
+
+  fSixGevData = 0; // Set 6 GeV data parameter to false unless set in parameter file
+  fDebugAdc   = 0; // Set ADC debug parameter to false unless set in parameter file
+
+  gHcParms->LoadParmValues((DBRequest*)&list, prefix);
+
+  if (fSixGevData) cout << "6 GeV Data Analysis Flag Set To TRUE" << endl;
 
   fIsInit = true;
 
+  cout << "Track Matching Parameters for: " << GetApparatus()->GetName()
+       << "." << GetName() << endl;
+  for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
+    cout << "Region = " << iregion + 1 << endl;
+    for (Int_t ivalue = 0; ivalue < 8; ivalue++)
+      cout << fRegionValue[GetIndex(iregion, ivalue)] << "  ";
+    cout << endl;
+  }
+
   return kOK;
 }
 
@@ -280,7 +340,7 @@ Int_t THcAerogel::DefineVariables( EMode mode )
 {
   // Initialize global variables for histogramming and tree
 
-  cout << "THcAerogel::DefineVariables called " << GetName() << endl;
+  cout << "THcAerogel::DefineVariables called for: " << GetName() << endl;
 
   if( mode == kDefine && fIsSetup ) return kOK;
   fIsSetup = ( mode == kDefine );
@@ -290,100 +350,119 @@ Int_t THcAerogel::DefineVariables( EMode mode )
   // Do we need to put the number of pos/neg TDC/ADC hits into the variables?
   // No.  They show up in tree as Ndata.H.aero.postdchits for example
 
-  RVarDef vars[] = {
-    {"postdchits", "List of Positive TDC hits",
-     "fPosTDCHits.THcSignalHit.GetPaddleNumber()"},
-    {"negtdchits", "List of Negative TDC hits",
-     "fNegTDCHits.THcSignalHit.GetPaddleNumber()"},
-    {"posadchits", "List of Positive ADC hits",
-     "fPosADCHits.THcSignalHit.GetPaddleNumber()"},
-    {"negadchits", "List of Negative ADC hits",
-     "fNegADCHits.THcSignalHit.GetPaddleNumber()"},
-    {"apos",  "Raw Positive ADC Amplitudes",   "fA_Pos"},
-    {"aneg",  "Raw Negative ADC Amplitudes",   "fA_Neg"},
-    {"apos_p",  "Ped-subtracted Positive ADC Amplitudes",   "fA_Pos_p"},
-    {"aneg_p",  "Ped-subtracted Negative ADC Amplitudes",   "fA_Neg_p"},
-    {"tpos",  "Raw Positive TDC",   "fT_Pos"},
-    {"tneg",  "Raw Negative TDC",   "fT_Neg"},
-    {"pos_npe","PEs Positive Tube","fPosNpe"},
-    {"neg_npe","PEs Negative Tube","fNegNpe"},
-    {"pos_npe_sum", "Total Positive Tube PEs", "fPosNpeSum"},
-    {"neg_npe_sum", "Total Negative Tube PEs", "fNegNpeSum"},
-    {"npe_sum", "Total PEs", "fNpeSum"},
-    {"ntdc_pos_hits", "Number of Positive Tube Hits", "fNTDCPosHits"},
-    {"ntdc_neg_hits", "Number of Negative Tube Hits", "fNTDCNegHits"},
-    {"ngood_hits", "Total number of good hits", "fNGoodHits"},
-
-    {"posGain", "List of positive PMT gains.", "fPosGain"},
-    {"negGain", "List of negative PMT gains.", "fNegGain"},
-
-    {"posAdcCounter",      "List of positive ADC counter numbers.",      "frPosAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
-    {"negAdcCounter",      "List of negative ADC counter numbers.",      "frNegAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
-
-    {"posAdcPedRaw",       "List of positive raw ADC pedestals",         "frPosAdcPedRaw.THcSignalHit.GetData()"},
-    {"posAdcPulseIntRaw",  "List of positive raw ADC pulse integrals.",  "frPosAdcPulseIntRaw.THcSignalHit.GetData()"},
-    {"posAdcPulseAmpRaw",  "List of positive raw ADC pulse amplitudes.", "frPosAdcPulseAmpRaw.THcSignalHit.GetData()"},
-    {"posAdcPulseTimeRaw", "List of positive raw ADC pulse times.",      "frPosAdcPulseTimeRaw.THcSignalHit.GetData()"},
-
-    {"posAdcPed",          "List of positive ADC pedestals",             "frPosAdcPed.THcSignalHit.GetData()"},
-    {"posAdcPulseInt",     "List of positive ADC pulse integrals.",      "frPosAdcPulseInt.THcSignalHit.GetData()"},
-    {"posAdcPulseAmp",     "List of positive ADC pulse amplitudes.",     "frPosAdcPulseAmp.THcSignalHit.GetData()"},
-
-    {"negAdcPedRaw",       "List of negative raw ADC pedestals",         "frNegAdcPedRaw.THcSignalHit.GetData()"},
-    {"negAdcPulseIntRaw",  "List of negative raw ADC pulse integrals.",  "frNegAdcPulseIntRaw.THcSignalHit.GetData()"},
-    {"negAdcPulseAmpRaw",  "List of negative raw ADC pulse amplitudes.", "frNegAdcPulseAmpRaw.THcSignalHit.GetData()"},
-    {"negAdcPulseTimeRaw", "List of negative raw ADC pulse times.",      "frNegAdcPulseTimeRaw.THcSignalHit.GetData()"},
-
-    {"negAdcPed",          "List of negative ADC pedestals",             "frNegAdcPed.THcSignalHit.GetData()"},
-    {"negAdcPulseInt",     "List of negative ADC pulse integrals.",      "frNegAdcPulseInt.THcSignalHit.GetData()"},
-    {"negAdcPulseAmp",     "List of negative ADC pulse amplitudes.",     "frNegAdcPulseAmp.THcSignalHit.GetData()"},
-
-    { 0 }
-  };
+  vector<RVarDef> vars;
+
+  vars.push_back({"posAdcCounter",   "Positive ADC counter numbers",   "frPosAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"});
+  vars.push_back({"negAdcCounter",   "Negative ADC counter numbers",   "frNegAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"});
+  vars.push_back({"posAdcErrorFlag", "Error Flag for When FPGA Fails", "fPosAdcErrorFlag.THcSignalHit.GetData()"});
+  vars.push_back({"negAdcErrorFlag", "Error Flag for When FPGA Fails", "fNegAdcErrorFlag.THcSignalHit.GetData()"});
+
+  vars.push_back({"numGoodPosAdcHits",    "Number of Good Positive ADC Hits Per PMT", "fNumGoodPosAdcHits"});    // Aerogel occupancy
+  vars.push_back({"numGoodNegAdcHits",    "Number of Good Negative ADC Hits Per PMT", "fNumGoodNegAdcHits"});    // Aerogel occupancy
+  vars.push_back({"totNumGoodPosAdcHits", "Total Number of Good Positive ADC Hits",   "fTotNumGoodPosAdcHits"}); // Aerogel multiplicity
+  vars.push_back({"totNumGoodNegAdcHits", "Total Number of Good Negative ADC Hits",   "fTotNumGoodNegAdcHits"}); // Aerogel multiplicity
+
+  vars.push_back({"totnumGoodAdcHits",   "TotalNumber of Good ADC Hits Per PMT",      "fTotNumGoodAdcHits"});    // Aerogel multiplicity
+  vars.push_back({"numTracksMatched",    "Number of Tracks Matched Per Region",       "fNumTracksMatched"});
+  vars.push_back({"numTracksFired",      "Number of Tracks that Fired Per Region",    "fNumTracksFired"});
+  vars.push_back({"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"});
+  vars.push_back({"totNumTracksFired",   "Total Number of Tracks that Fired",         "fTotNumTracksFired"});
+
+  vars.push_back({"posNpe",    "Number of Positive PEs",       "fPosNpe"});
+  vars.push_back({"negNpe",    "Number of Negative PEs",       "fNegNpe"});
+  vars.push_back({"posNpeSum", "Total Number of Positive PEs", "fPosNpeSum"});
+  vars.push_back({"negNpeSum", "Total Number of Negative PEs", "fNegNpeSum"});
+  vars.push_back({"npeSum",    "Total Number of PEs",          "fNpeSum"});
+
+  vars.push_back({"goodPosAdcPed",         "Good Negative ADC pedestals",           "fGoodPosAdcPed"});
+  vars.push_back({"goodPosAdcPulseInt",    "Good Negative ADC pulse integrals",     "fGoodPosAdcPulseInt"});
+  vars.push_back({"goodPosAdcPulseIntRaw", "Good Negative ADC raw pulse integrals", "fGoodPosAdcPulseIntRaw"});
+  vars.push_back({"goodPosAdcPulseAmp",    "Good Negative ADC pulse amplitudes",    "fGoodPosAdcPulseAmp"});
+  vars.push_back({"goodPosAdcPulseTime",   "Good Negative ADC pulse times",         "fGoodPosAdcPulseTime"});
+
+  vars.push_back({"goodNegAdcPed",         "Good Negative ADC pedestals",           "fGoodNegAdcPed"});
+  vars.push_back({"goodNegAdcPulseInt",    "Good Negative ADC pulse integrals",     "fGoodNegAdcPulseInt"});
+  vars.push_back({"goodNegAdcPulseIntRaw", "Good Negative ADC raw pulse integrals", "fGoodNegAdcPulseIntRaw"});
+  vars.push_back({"goodNegAdcPulseAmp",    "Good Negative ADC pulse amplitudes",    "fGoodNegAdcPulseAmp"});
+  vars.push_back({"goodNegAdcPulseTime",   "Good Negative ADC pulse times",         "fGoodNegAdcPulseTime"});
+
+  if (fDebugAdc) {
+    vars.push_back({"posGain", "Positive PMT gains", "fPosGain"});
+    vars.push_back({"negGain", "Negative PMT gains", "fNegGain"});
+
+    vars.push_back({"numPosAdcHits",        "Number of Positive ADC Hits Per PMT",      "fNumPosAdcHits"});        // Aerogel occupancy
+    vars.push_back({"totNumPosAdcHits",     "Total Number of Positive ADC Hits",        "fTotNumPosAdcHits"});     // Aerogel multiplicity
+    vars.push_back({"numNegAdcHits",        "Number of Negative ADC Hits Per PMT",      "fNumNegAdcHits"});        // Aerogel occupancy
+    vars.push_back({"totNumNegAdcHits",     "Total Number of Negative ADC Hits",        "fTotNumNegAdcHits"});     // Aerogel multiplicity
+    vars.push_back({"totnumAdcHits",       "Total Number of ADC Hits Per PMT",          "fTotNumAdcHits"});        // Aerogel multiplicity
+
+    vars.push_back({"posAdcPedRaw",       "Positive Raw ADC pedestals",        "frPosAdcPedRaw.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPulseIntRaw",  "Positive Raw ADC pulse integrals",  "frPosAdcPulseIntRaw.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPulseAmpRaw",  "Positive Raw ADC pulse amplitudes", "frPosAdcPulseAmpRaw.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPulseTimeRaw", "Positive Raw ADC pulse times",      "frPosAdcPulseTimeRaw.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPed",          "Positive ADC pedestals",            "frPosAdcPed.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPulseInt",     "Positive ADC pulse integrals",      "frPosAdcPulseInt.THcSignalHit.GetData()"});
+    vars.push_back({"posAdcPulseAmp",     "Positive ADC pulse amplitudes",     "frPosAdcPulseAmp.THcSignalHit.GetData()"});
+
+    vars.push_back({"negAdcPedRaw",       "Negative Raw ADC pedestals",        "frNegAdcPedRaw.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPulseIntRaw",  "Negative Raw ADC pulse integrals",  "frNegAdcPulseIntRaw.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPulseAmpRaw",  "Negative Raw ADC pulse amplitudes", "frNegAdcPulseAmpRaw.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPulseTimeRaw", "Negative Raw ADC pulse times",      "frNegAdcPulseTimeRaw.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPed",          "Negative ADC pedestals",            "frNegAdcPed.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPulseInt",     "Negative ADC pulse integrals",      "frNegAdcPulseInt.THcSignalHit.GetData()"});
+    vars.push_back({"negAdcPulseAmp",     "Negative ADC pulse amplitudes",     "frNegAdcPulseAmp.THcSignalHit.GetData()"});
+  }
+
+  if (fSixGevData) {
+    vars.push_back({"apos",            "Positive Raw ADC Amplitudes",            "fA_Pos"});
+    vars.push_back({"aneg",            "Negative Raw ADC Amplitudes",            "fA_Neg"});
+    vars.push_back({"apos_p",          "Positive Ped-subtracted ADC Amplitudes", "fA_Pos_p"});
+    vars.push_back({"aneg_p",          "Negative Ped-subtracted ADC Amplitudes", "fA_Neg_p"});
+    vars.push_back({"tpos",            "Positive Raw TDC",                       "fT_Pos"});
+    vars.push_back({"tneg",            "Negative Raw TDC",                       "fT_Neg"});
+    vars.push_back({"ntdc_pos_hits",   "Number of Positive Tube Hits",           "fNTDCPosHits"});
+    vars.push_back({"ntdc_neg_hits",   "Number of Negative Tube Hits",           "fNTDCNegHits"});
+    vars.push_back({"posadchits",      "Positive ADC hits",                      "fPosADCHits.THcSignalHit.GetPaddleNumber()"});
+    vars.push_back({"negadchits",      "Negative ADC hits",                      "fNegADCHits.THcSignalHit.GetPaddleNumber()"});
+    vars.push_back({"postdchits",      "Positive TDC hits",                      "fPosTDCHits.THcSignalHit.GetPaddleNumber()"});
+    vars.push_back({"negtdchits",      "Negative TDC hits",                      "fNegTDCHits.THcSignalHit.GetPaddleNumber()"});
+    vars.push_back({"nGoodHits",       "Total number of good hits",              "fNGoodHits"});
+    vars.push_back({"posNpeSixGev",    "Number of Positive PEs",                 "fPosNpeSixGev"});
+    vars.push_back({"negNpeSixGev",    "Number of Negative PEs",                 "fNegNpeSixGev"});
+    vars.push_back({"posNpeSumSixGev", "Total Number of Positive PEs",           "fPosNpeSumSixGev"});
+    vars.push_back({"negNpeSumSixGev", "Total Number of Negative PEs",           "fNegNpeSumSixGev"});
+    vars.push_back({"npeSumSixGev",    "Total Number of PEs",                    "fNpeSumSixGev"});
+  }
+
+  RVarDef end {0};
+  vars.push_back(end);
+
+  return DefineVarsFromList(vars.data(), mode);
 
-  return DefineVarsFromList( vars, mode );
 }
 //_____________________________________________________________________________
 inline
 void THcAerogel::Clear(Option_t* opt)
 {
   // Clear the hit lists
-  fPosTDCHits->Clear();
-  fNegTDCHits->Clear();
-  fPosADCHits->Clear();
-  fNegADCHits->Clear();
-
-  // Clear Aerogel variables  from h_aero.f
-
-  fNhits = 0;	     // Don't really need to do this.  (Be sure this called before Decode)
-
+  fNhits                = 0;
+  fTotNumAdcHits        = 0;
+  fTotNumGoodAdcHits    = 0;
+  fTotNumPosAdcHits     = 0;
+  fTotNumNegAdcHits     = 0;
+  fTotNumGoodPosAdcHits = 0;
+  fTotNumGoodNegAdcHits = 0;
+  fTotNumTracksMatched  = 0;
+  fTotNumTracksFired    = 0;
+
+  fNpeSum    = 0.0;
   fPosNpeSum = 0.0;
   fNegNpeSum = 0.0;
-  fNpeSum = 0.0;
-
-  fNGoodHits = 0;
-
-  fNADCPosHits = 0;
-  fNADCNegHits = 0;
-  fNTDCPosHits = 0;
-  fNTDCNegHits = 0;
-
-  for(Int_t itube = 0;itube < fNelem;itube++) {
-    fA_Pos[itube] = 0;
-    fA_Neg[itube] = 0;
-    fA_Pos_p[itube] = 0;
-    fA_Neg_p[itube] = 0;
-    fT_Pos[itube] = 0;
-    fT_Neg[itube] = 0;
-    fPosNpe[itube] = 0.0;
-    fNegNpe[itube] = 0.0;
-  }
 
   frPosAdcPedRaw->Clear();
   frPosAdcPulseIntRaw->Clear();
   frPosAdcPulseAmpRaw->Clear();
   frPosAdcPulseTimeRaw->Clear();
-
   frPosAdcPed->Clear();
   frPosAdcPulseInt->Clear();
   frPosAdcPulseAmp->Clear();
@@ -392,209 +471,146 @@ void THcAerogel::Clear(Option_t* opt)
   frNegAdcPulseIntRaw->Clear();
   frNegAdcPulseAmpRaw->Clear();
   frNegAdcPulseTimeRaw->Clear();
-
   frNegAdcPed->Clear();
   frNegAdcPulseInt->Clear();
   frNegAdcPulseAmp->Clear();
-}
 
-//_____________________________________________________________________________
-Int_t THcAerogel::Decode( const THaEvData& evdata )
-{
-  // Get the Hall C style hitlist (fRawHitList) for this event
-  fNhits = DecodeToHitList(evdata);
+  fPosAdcErrorFlag->Clear();
+  fNegAdcErrorFlag->Clear();
+
+  for (UInt_t ielem = 0; ielem < fNumPosAdcHits.size(); ielem++)
+    fNumPosAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumNegAdcHits.size(); ielem++)
+    fNumNegAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumGoodPosAdcHits.size(); ielem++)
+    fNumGoodPosAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumGoodNegAdcHits.size(); ielem++)
+    fNumGoodNegAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumTracksMatched.size(); ielem++)
+    fNumTracksMatched.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumTracksFired.size(); ielem++)
+    fNumTracksFired.at(ielem) = 0;
+
+  for (UInt_t ielem = 0; ielem < fGoodPosAdcPed.size(); ielem++) {
+    fGoodPosAdcPed.at(ielem)         = 0.0;
+    fGoodPosAdcPulseInt.at(ielem)    = 0.0;
+    fGoodPosAdcPulseIntRaw.at(ielem) = 0.0;
+    fGoodPosAdcPulseAmp.at(ielem)    = 0.0;
+    fGoodPosAdcPulseTime.at(ielem)   = 0.0;
+    fPosNpe.at(ielem)                = 0.0;
+  }
+  for (UInt_t ielem = 0; ielem < fGoodNegAdcPed.size(); ielem++) {
+    fGoodNegAdcPed.at(ielem)         = 0.0;
+    fGoodNegAdcPulseInt.at(ielem)    = 0.0;
+    fGoodNegAdcPulseIntRaw.at(ielem) = 0.0;
+    fGoodNegAdcPulseAmp.at(ielem)    = 0.0;
+    fGoodNegAdcPulseTime.at(ielem)   = 0.0;
+    fNegNpe.at(ielem)                = 0.0;
+  }
 
-  if(gHaCuts->Result("Pedestal_event")) {
+  // 6 GeV variables
+  fNGoodHits       = 0;
+  fNADCPosHits     = 0;
+  fNADCNegHits     = 0;
+  fNTDCPosHits     = 0;
+  fNTDCNegHits     = 0;
+  fNpeSumSixGev    = 0.0;
+  fPosNpeSumSixGev = 0.0;
+  fNegNpeSumSixGev = 0.0;
+  fPosTDCHits->Clear();
+  fNegTDCHits->Clear();
+  fPosADCHits->Clear();
+  fNegADCHits->Clear();
 
-    AccumulatePedestals(fRawHitList);
+  for (UInt_t ielem = 0; ielem < fPosNpeSixGev.size(); ielem++)
+    fPosNpeSixGev.at(ielem) = 0.0;
+  for (UInt_t ielem = 0; ielem < fNegNpeSixGev.size(); ielem++)
+    fNegNpeSixGev.at(ielem) = 0.0;
 
-    fAnalyzePedestals = 1;	// Analyze pedestals first normal events
-    return(0);
+  for(Int_t itube = 0;itube < fNelem;itube++) {
+    fA_Pos[itube]   = 0;
+    fA_Neg[itube]   = 0;
+    fA_Pos_p[itube] = 0;
+    fA_Neg_p[itube] = 0;
+    fT_Pos[itube]   = 0;
+    fT_Neg[itube]   = 0;
   }
 
-  if(fAnalyzePedestals) {
+}
 
-    CalculatePedestals();
-    Print("");
+//_____________________________________________________________________________
+Int_t THcAerogel::Decode( const THaEvData& evdata )
+{
+  // Get the Hall C style hitlist (fRawHitList) for this event
+  fNhits = DecodeToHitList(evdata);
 
-    fAnalyzePedestals = 0;	// Don't analyze pedestals next event
+  if (fSixGevData) {
+    if(gHaCuts->Result("Pedestal_event")) {
+      AccumulatePedestals(fRawHitList);
+      fAnalyzePedestals = 1;	// Analyze pedestals first normal events
+      return(0);
+    }
+    if(fAnalyzePedestals) {
+      CalculatePedestals();
+      Print("");
+      fAnalyzePedestals = 0;	// Don't analyze pedestals next event
+    }
   }
 
-
-  Int_t ihit = 0;
-  Int_t nPosTDCHits=0;
-  Int_t nNegTDCHits=0;
-  Int_t nPosADCHits=0;
-  Int_t nNegADCHits=0;
-
+  Int_t  ihit         = 0;
   UInt_t nrPosAdcHits = 0;
   UInt_t nrNegAdcHits = 0;
 
   while(ihit < fNhits) {
-    THcAerogelHit* hit = (THcAerogelHit *) fRawHitList->At(ihit);
-
-    Int_t padnum = hit->fCounter;
+    THcAerogelHit* hit          = (THcAerogelHit*) fRawHitList->At(ihit);
+    Int_t          npmt         = hit->fCounter;
+    THcRawAdcHit&  rawPosAdcHit = hit->GetRawAdcHitPos();
+    THcRawAdcHit&  rawNegAdcHit = hit->GetRawAdcHitNeg();
 
-    THcRawAdcHit& rawPosAdcHit = hit->GetRawAdcHitPos();
     for (UInt_t thit=0; thit<rawPosAdcHit.GetNPulses(); ++thit) {
-      ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPedRaw());
-      ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPed());
-
-      ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseIntRaw(thit));
-      ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseInt(thit));
-
-      ((THcSignalHit*) frPosAdcPulseAmpRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseAmpRaw(thit));
-      ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseAmp(thit));
-
-      ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(nrPosAdcHits))->Set(padnum, rawPosAdcHit.GetPulseTimeRaw(thit));
-
-      ++nrPosAdcHits;
-    }
-    THcRawAdcHit& rawNegAdcHit = hit->GetRawAdcHitNeg();
-    for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
-      ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPedRaw());
-      ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPed());
-
-      ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseIntRaw(thit));
-      ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseInt(thit));
-
-      ((THcSignalHit*) frNegAdcPulseAmpRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseAmpRaw(thit));
-      ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseAmp(thit));
 
-      ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(nrNegAdcHits))->Set(padnum, rawNegAdcHit.GetPulseTimeRaw(thit));
+      ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPedRaw());
+      ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPed());
 
-      ++nrNegAdcHits;
-    }
+      ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseIntRaw(thit));
+      ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseInt(thit));
 
-    Int_t adc_pos;
-    Int_t adc_neg;
-    Int_t tdc_pos=-1;
-    Int_t tdc_neg=-1;
-   // TDC positive hit
-    if(hit->GetRawTdcHitPos().GetNHits() >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
-      tdc_pos = hit->GetRawTdcHitPos().GetTime()+fTdcOffset;
-      sighit->Set(hit->fCounter, tdc_pos);
-    }
-
-    // TDC negative hit
-    if(hit->GetRawTdcHitNeg().GetNHits() >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
-      tdc_neg = hit->GetRawTdcHitNeg().GetTime()+fTdcOffset;
-      sighit->Set(hit->fCounter, tdc_neg);
-    }
+      ((THcSignalHit*) frPosAdcPulseAmpRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseAmpRaw(thit));
+      ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseAmp(thit));
 
-    // ADC positive hit
-    if((adc_pos = hit->GetRawAdcHitPos().GetPulseInt()) > 0) {
-      THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
-      sighit->Set(hit->fCounter, adc_pos);
-    }
+      ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, rawPosAdcHit.GetPulseTimeRaw(thit));
 
-    // ADC negative hit
-    if((adc_neg = hit->GetRawAdcHitNeg().GetPulseInt()) > 0) {
-      THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
-      sighit->Set(hit->fCounter, adc_neg);
-    }
+      if (rawPosAdcHit.GetPulseAmpRaw(thit) > 0)  ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 0);
+      if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 1);
 
-    // For each TDC, identify the first hit that is positive.
-    tdc_pos = -1;
-    tdc_neg = -1;
-    for(UInt_t thit=0; thit<hit->GetRawTdcHitPos().GetNHits(); thit++) {
-      Int_t tdc = hit->GetRawTdcHitPos().GetTime(thit);
-      if(tdc >=0 ) {
-	tdc_pos = tdc;
-	break;
-      }
-    }
-    for(UInt_t thit=0; thit<hit->GetRawTdcHitNeg().GetNHits(); thit++) {
-      Int_t tdc = hit->GetRawTdcHitNeg().GetTime(thit);
-      if(tdc >= 0) {
-	tdc_neg = tdc;
-	break;
-      }
+      ++nrPosAdcHits;
+      fTotNumAdcHits++;
+      fTotNumPosAdcHits++;
+      fNumPosAdcHits.at(npmt-1) = npmt;
     }
 
-    // Fill the the per detector ADC and TDC arrays
-    Int_t npmt = hit->fCounter - 1;
+    for (UInt_t thit=0; thit<rawNegAdcHit.GetNPulses(); ++thit) {
+      ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPedRaw());
+      ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPed());
 
-    fA_Pos[npmt] = adc_pos;
-    fA_Neg[npmt] = adc_neg;
-    fA_Pos_p[npmt] = fA_Pos[npmt] - fPosPedMean[npmt];
-    fA_Neg_p[npmt] = fA_Neg[npmt] - fNegPedMean[npmt];
-    fT_Pos[npmt] = tdc_pos;
-    fT_Neg[npmt] = tdc_neg;
+      ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseIntRaw(thit));
+      ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseInt(thit));
 
-    if(fA_Pos[npmt] < 8000) {
-      fPosNpe[npmt] = fPosGain[npmt]*fA_Pos_p[npmt];
-    } else {
-      fPosNpe[npmt] = 100.0;
-    }
+      ((THcSignalHit*) frNegAdcPulseAmpRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseAmpRaw(thit));
+      ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseAmp(thit));
 
-    if(fA_Neg[npmt] < 8000) {
-      fNegNpe[npmt] = fNegGain[npmt]*fA_Neg_p[npmt];
-    } else {
-      fNegNpe[npmt] = 100.0;
-    }
+      ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, rawNegAdcHit.GetPulseTimeRaw(thit));
 
-    fPosNpeSum += fPosNpe[npmt];
-    fNegNpeSum += fNegNpe[npmt];
+      if (rawNegAdcHit.GetPulseAmpRaw(thit) > 0)  ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 0);
+      if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 1);
 
-    // Sum positive and negative hits to fill tot_good_hits
-    if(fPosNpe[npmt] > 0.3) {
-      fNADCPosHits++;
-      fNGoodHits++;
-    }
-    if(fNegNpe[npmt] > 0.3) {
-      fNADCNegHits++;
-      fNGoodHits++;
-    }
-    if(fT_Pos[npmt] > 0 && fT_Pos[npmt] < 8000) {
-      fNTDCPosHits++;
-    }
-    if(fT_Neg[npmt] > 0 && fT_Neg[npmt] < 8000) {
-      fNTDCNegHits++;
+      ++nrNegAdcHits;
+      fTotNumAdcHits++;
+      fTotNumNegAdcHits++;
+      fNumNegAdcHits.at(npmt-1) = npmt;
     }
-
     ihit++;
   }
-
-  if(fPosNpeSum > 0.5 || fNegNpeSum > 0.5) {
-    fNpeSum = fPosNpeSum + fNegNpeSum;
-  } else {
-    fNpeSum = 0.0;
-  }
-
-  // If total hits are 0, then give a noticable ridiculous NPE
-  if(fNhits < 1) {
-    fNpeSum = 0.0;
-  }
-
-  // The following code is in the fortran.  It probably doesn't work
-  // right because the arrays are not cleared first and the aero_ep,
-  // aero_en, ... lines make no sense.
-
-  //* Next, fill the rawadc variables with the actual tube values
-  //*       mainly for diagnostic purposes.
-  //
-  //      do ihit=1,haero_tot_hits
-  //
-  //         npmt=haero_pair_num(ihit)
-  //
-  //         haero_rawadc_pos(npmt)=haero_adc_pos(ihit)
-  //         aero_ep(npmt)=haero_rawadc_pos(ihit)
-  //
-  //         haero_rawadc_neg(npmt)=haero_adc_neg(ihit)
-  //         aero_en(npmt)=haero_rawadc_neg(ihit)
-  //
-  //         haero_rawtdc_neg(npmt)=haero_tdc_neg(ihit)
-  //         aero_tn(npmt)= haero_tdc_neg(ihit)
-  //
-  //         haero_rawtdc_pos(npmt)=haero_tdc_pos(ihit)
-  //         aero_tp(npmt)= haero_tdc_pos(ihit)
-  //
-  //      enddo
-
-
   return ihit;
 }
 
@@ -608,96 +624,300 @@ Int_t THcAerogel::ApplyCorrections( void )
 Int_t THcAerogel::CoarseProcess( TClonesArray&  ) //tracks
 {
 
-  // All code previously here moved into decode
+    // Loop over the elements in the TClonesArray
+    for(Int_t ielem = 0; ielem < frPosAdcPulseInt->GetEntries(); ielem++) {
+
+      Int_t    npmt         = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+      Double_t pulsePed     = ((THcSignalHit*) frPosAdcPed->ConstructedAt(ielem))->GetData();
+      Double_t pulseInt     = ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(ielem))->GetData();
+      Double_t pulseIntRaw  = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+      Double_t pulseAmp     = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(ielem))->GetData();
+      Double_t pulseTime    = ((THcSignalHit*) frPosAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+      Bool_t   errorFlag    = ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
+      Bool_t   pulseTimeCut = pulseTime > fAdcTimeWindowMin && pulseTime < fAdcTimeWindowMax;
+
+      // By default, the last hit within the timing cut will be considered "good"
+      if (!errorFlag && pulseTimeCut) {
+    	fGoodPosAdcPed.at(npmt)         = pulsePed;
+    	fGoodPosAdcPulseInt.at(npmt)    = pulseInt;
+    	fGoodPosAdcPulseIntRaw.at(npmt) = pulseIntRaw;
+    	fGoodPosAdcPulseAmp.at(npmt)    = pulseAmp;
+    	fGoodPosAdcPulseTime.at(npmt)   = pulseTime;
+
+    	fPosNpe.at(npmt) = fPosGain[npmt]*fGoodPosAdcPulseInt.at(npmt);
+	fPosNpeSum += fPosNpe.at(npmt);
+    	fNpeSum += fPosNpeSum;
+
+	fTotNumGoodAdcHits++;
+    	fTotNumGoodPosAdcHits++;
+    	fNumGoodPosAdcHits.at(npmt) = npmt + 1;
+      }
+    }
 
-  ApplyCorrections();
+    // Loop over the elements in the TClonesArray
+    for(Int_t ielem = 0; ielem < frNegAdcPulseInt->GetEntries(); ielem++) {
+
+      Int_t    npmt         = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+      Double_t pulsePed     = ((THcSignalHit*) frNegAdcPed->ConstructedAt(ielem))->GetData();
+      Double_t pulseInt     = ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(ielem))->GetData();
+      Double_t pulseIntRaw  = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+      Double_t pulseAmp     = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
+      Double_t pulseTime    = ((THcSignalHit*) frNegAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+      Bool_t   errorFlag    = ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
+      Bool_t   pulseTimeCut = pulseTime > fAdcTimeWindowMin && pulseTime < fAdcTimeWindowMax;
+
+      // By default, the last hit within the timing cut will be considered "good"
+      if (!errorFlag && pulseTimeCut) {
+    	fGoodNegAdcPed.at(npmt)         = pulsePed;
+    	fGoodNegAdcPulseInt.at(npmt)    = pulseInt;
+    	fGoodNegAdcPulseIntRaw.at(npmt) = pulseIntRaw;
+    	fGoodNegAdcPulseAmp.at(npmt)    = pulseAmp;
+    	fGoodNegAdcPulseTime.at(npmt)   = pulseTime;
+
+    	fNegNpe.at(npmt) = fNegGain[npmt]*fGoodNegAdcPulseInt.at(npmt);
+	fNegNpeSum += fNegNpe.at(npmt);
+    	fNpeSum += fNegNpeSum;
+
+	fTotNumGoodAdcHits++;
+   	fTotNumGoodNegAdcHits++;
+    	fNumGoodNegAdcHits.at(npmt) = npmt + 1;
+      }
+    }
 
-  return 0;
+    for(Int_t ihit=0; ihit < fNhits; ihit++) {
+
+      Int_t nPosTDCHits = 0;
+      Int_t nNegTDCHits = 0;
+      Int_t nPosADCHits = 0;
+      Int_t nNegADCHits = 0;
+
+      // 6 GeV calculations
+      Int_t adc_pos;
+      Int_t adc_neg;
+      Int_t tdc_pos = -1;
+      Int_t tdc_neg = -1;
+      if (fSixGevData) {
+	THcAerogelHit* hit = (THcAerogelHit*) fRawHitList->At(ihit);
+	Int_t npmt = hit->fCounter - 1;
+
+	// Sum positive and negative hits to fill tot_good_hits
+	if(fPosNpe.at(npmt) > 0.3) {fNADCPosHits++; fNGoodHits++;}
+	if(fNegNpe.at(npmt) > 0.3) {fNADCNegHits++; fNGoodHits++;}
+
+	// ADC positive hit
+	if((adc_pos = hit->GetRawAdcHitPos().GetPulseInt()) > 0) {
+	  THcSignalHit *sighit = (THcSignalHit*) fPosADCHits->ConstructedAt(nPosADCHits++);
+	  sighit->Set(hit->fCounter, adc_pos);
+	}
+	// ADC negative hit
+	if((adc_neg = hit->GetRawAdcHitNeg().GetPulseInt()) > 0) {
+	  THcSignalHit *sighit = (THcSignalHit*) fNegADCHits->ConstructedAt(nNegADCHits++);
+	  sighit->Set(hit->fCounter, adc_neg);
+	}
+	// TDC positive hit
+	if(hit->GetRawTdcHitPos().GetNHits() >  0) {
+	  THcSignalHit *sighit = (THcSignalHit*) fPosTDCHits->ConstructedAt(nPosTDCHits++);
+	  tdc_pos = hit->GetRawTdcHitPos().GetTime()+fTdcOffset;
+	  sighit->Set(hit->fCounter, tdc_pos);
+	}
+	// TDC negative hit
+	if(hit->GetRawTdcHitNeg().GetNHits() >  0) {
+	  THcSignalHit *sighit = (THcSignalHit*) fNegTDCHits->ConstructedAt(nNegTDCHits++);
+	  tdc_neg = hit->GetRawTdcHitNeg().GetTime()+fTdcOffset;
+	  sighit->Set(hit->fCounter, tdc_neg);
+	}
+	// For each TDC, identify the first hit that is positive.
+	tdc_pos = -1;
+	tdc_neg = -1;
+	for(UInt_t thit=0; thit<hit->GetRawTdcHitPos().GetNHits(); thit++) {
+	  Int_t tdc = hit->GetRawTdcHitPos().GetTime(thit);
+	  if(tdc >=0 ) {
+	    tdc_pos = tdc;
+	    break;
+	  }
+	}
+	for(UInt_t thit=0; thit<hit->GetRawTdcHitNeg().GetNHits(); thit++) {
+	  Int_t tdc = hit->GetRawTdcHitNeg().GetTime(thit);
+	  if(tdc >= 0) {
+	    tdc_neg = tdc;
+	    break;
+	  }
+	}
+
+	fA_Pos[npmt] = adc_pos;
+	fA_Neg[npmt] = adc_neg;
+	fA_Pos_p[npmt] = fA_Pos[npmt] - fPosPedMean[npmt];
+	fA_Neg_p[npmt] = fA_Neg[npmt] - fNegPedMean[npmt];
+	fT_Pos[npmt] = tdc_pos;
+	fT_Neg[npmt] = tdc_neg;
+
+	if(fA_Pos[npmt] < 8000) fPosNpeSixGev[npmt] = fPosGain[npmt]*fA_Pos_p[npmt];
+	else fPosNpeSixGev[npmt] = 100.0;
+
+	if(fA_Neg[npmt] < 8000) fNegNpeSixGev[npmt] = fNegGain[npmt]*fA_Neg_p[npmt];
+	else fNegNpeSixGev[npmt] = 100.0;
+
+	fPosNpeSumSixGev += fPosNpeSixGev[npmt];
+	fNegNpeSumSixGev += fNegNpeSixGev[npmt];
+
+	// Sum positive and negative hits to fill tot_good_hits
+	if(fPosNpeSixGev[npmt] > 0.3) {fNADCPosHits++; fNGoodHits++;}
+	if(fNegNpeSixGev[npmt] > 0.3) {fNADCNegHits++; fNGoodHits++;}
+
+	if(fT_Pos[npmt] > 0 && fT_Pos[npmt] < 8000) fNTDCPosHits++;
+	if(fT_Neg[npmt] > 0 && fT_Neg[npmt] < 8000) fNTDCNegHits++;
+      }
+
+      if (fPosNpeSumSixGev > 0.5 || fNegNpeSumSixGev > 0.5)
+	fNpeSumSixGev = fPosNpeSumSixGev + fNegNpeSumSixGev;
+      else fNpeSumSixGev = 0.0;
+      // If total hits are 0, then give a noticable ridiculous NPE
+      if (fNhits < 1) fNpeSumSixGev = 0.0;
+
+    }
+    return 0;
 }
 
 //_____________________________________________________________________________
 Int_t THcAerogel::FineProcess( TClonesArray& tracks )
 {
 
+  Int_t nTracks = tracks.GetLast() + 1;
+
+  for (Int_t itrack = 0; itrack < nTracks; itrack++) {
+
+    THaTrack* track = dynamic_cast<THaTrack*> (tracks[itrack]);
+    if (track->GetIndex() != 0) continue;  // Select the best track
+
+    Double_t trackChi2    = track->GetChi2();
+    Int_t    trackNDoF    = track->GetNDoF();
+    Double_t trackRedChi2 = trackChi2/trackNDoF;
+    Double_t trackBeta    = track->GetBeta();
+    Double_t trackEnergy  = track->GetEnergy();
+    Double_t trackMom     = track->GetP();
+    Double_t trackENorm   = trackEnergy/trackMom;
+    Double_t trackXfp     = track->GetX();
+    Double_t trackYfp     = track->GetY();
+    Double_t trackTheta   = track->GetTheta();
+    Double_t trackPhi     = track->GetPhi();
+
+    Bool_t trackRedChi2Cut = trackRedChi2 > fRedChi2Min && trackRedChi2 < fRedChi2Max;
+    Bool_t trackBetaCut    = trackBeta    > fBetaMin    && trackBeta    < fBetaMax;
+    Bool_t trackENormCut   = trackENorm   > fENormMin   && trackENorm   < fENormMax;
+
+    if (trackRedChi2Cut && trackBetaCut && trackENormCut) {
+
+      // Project the track to the Aerogel diffuser box plane
+      Double_t xAtAero = trackXfp + trackTheta * fDiffBoxZPos;
+      Double_t yAtAero = trackYfp + trackPhi   * fDiffBoxZPos;
+
+      // cout << "Aerogel Detector: " << GetName() << endl;
+      // cout << "nTracks = " << nTracks << "\t" << "trackChi2 = " << trackChi2
+      // 	   << "\t" << "trackNDof = " << trackNDoF << "\t" << "trackRedChi2 = " << trackRedChi2 << endl;
+      // cout << "trackBeta = " << trackBeta << "\t" << "trackEnergy = " << trackEnergy << "\t"
+      // 	   << "trackMom = " << trackMom << "\t" << "trackENorm = " << trackENorm << endl;
+      // cout << "trackXfp = " << trackXfp << "\t" << "trackYfp = " << trackYfp << "\t"
+      // 	   << "trackTheta = " << trackTheta << "\t" << "trackPhi = " << trackPhi << endl;
+      // cout << "fDiffBoxZPos = " << fDiffBoxZPos << "\t" << "xAtAero = " << xAtAero << "\t" << "yAtAero = " << yAtAero << endl;
+      // cout << "=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:" << endl;
+
+
+      for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
+
+      	if ((TMath::Abs(fRegionValue[GetIndex(iregion, 0)] - xAtAero)    < fRegionValue[GetIndex(iregion, 4)]) &&
+      	    (TMath::Abs(fRegionValue[GetIndex(iregion, 1)] - yAtAero)    < fRegionValue[GetIndex(iregion, 5)]) &&
+      	    (TMath::Abs(fRegionValue[GetIndex(iregion, 2)] - trackTheta) < fRegionValue[GetIndex(iregion, 6)]) &&
+      	    (TMath::Abs(fRegionValue[GetIndex(iregion, 3)] - trackPhi)   < fRegionValue[GetIndex(iregion, 7)])) {
+
+	  fTotNumTracksMatched++;
+      	  fNumTracksMatched.at(iregion) = iregion + 1;
+
+      	  if (fNpeSum > fNpeThresh) {
+      	    fTotNumTracksFired++;
+      	    fNumTracksFired.at(iregion) = iregion + 1;
+      	  }  // NPE threshold cut
+      	}  // Regional cuts
+      }  // Loop over regions
+    }  // Tracking cuts
+  }  // Track loop
+
   return 0;
 }
 
 //_____________________________________________________________________________
-void THcAerogel::InitializePedestals( )
+// Method for initializing pedestals in the 6 GeV era
+void THcAerogel::InitializePedestals()
 {
   fNPedestalEvents = 0;
-  fMinPeds = 0;                    // Do not calculate pedestals by default
-  fPosPedSum = new Int_t [fNelem];
-  fPosPedSum2 = new Int_t [fNelem];
+  fMinPeds         = 0;                    // Do not calculate pedestals by default
+
+  fPosPedSum   = new Int_t [fNelem];
+  fPosPedSum2  = new Int_t [fNelem];
   fPosPedLimit = new Int_t [fNelem];
   fPosPedCount = new Int_t [fNelem];
-  fNegPedSum = new Int_t [fNelem];
-  fNegPedSum2 = new Int_t [fNelem];
+  fNegPedSum   = new Int_t [fNelem];
+  fNegPedSum2  = new Int_t [fNelem];
   fNegPedLimit = new Int_t [fNelem];
   fNegPedCount = new Int_t [fNelem];
-
-  fPosPed = new Double_t [fNelem];
-  fNegPed = new Double_t [fNelem];
-  fPosThresh = new Double_t [fNelem];
-  fNegThresh = new Double_t [fNelem];
-  for(Int_t i=0;i<fNelem;i++) {
-    fPosPedSum[i] = 0;
-    fPosPedSum2[i] = 0;
+  fPosPed      = new Double_t [fNelem];
+  fNegPed      = new Double_t [fNelem];
+  fPosThresh   = new Double_t [fNelem];
+  fNegThresh   = new Double_t [fNelem];
+
+  for(Int_t i = 0;i < fNelem; i++) {
+    fPosPedSum[i]   = 0;
+    fPosPedSum2[i]  = 0;
     fPosPedLimit[i] = 1000;   // In engine, this are set in parameter file
     fPosPedCount[i] = 0;
-    fNegPedSum[i] = 0;
-    fNegPedSum2[i] = 0;
+    fNegPedSum[i]   = 0;
+    fNegPedSum2[i]  = 0;
     fNegPedLimit[i] = 1000;   // In engine, this are set in parameter file
     fNegPedCount[i] = 0;
-    fPosPedMean[i] = 0;       // Default pedestal values
-    fNegPedMean[i] = 0;       // Default pedestal values
+    fPosPedMean[i]  = 0;       // Default pedestal values
+    fNegPedMean[i]  = 0;       // Default pedestal values
   }
 
-  fPosNpe = new Double_t [fNelem];
-  fNegNpe = new Double_t [fNelem];
 }
 
 //_____________________________________________________________________________
+// Method for accumulating pedestals in the 6 GeV era
 void THcAerogel::AccumulatePedestals(TClonesArray* rawhits)
 {
   // Extract data from the hit list, accumulating into arrays for
   // calculating pedestals
 
   Int_t nrawhits = rawhits->GetLast()+1;
+  Int_t ihit     = 0;
 
-  Int_t ihit = 0;
   while(ihit < nrawhits) {
     THcAerogelHit* hit = (THcAerogelHit *) rawhits->At(ihit);
 
     Int_t element = hit->fCounter - 1;
-    Int_t adcpos = hit->GetRawAdcHitPos().GetPulseInt();
-    Int_t adcneg = hit->GetRawAdcHitNeg().GetPulseInt();
+    Int_t adcpos  = hit->GetRawAdcHitPos().GetPulseInt();
+    Int_t adcneg  = hit->GetRawAdcHitNeg().GetPulseInt();
     if(adcpos <= fPosPedLimit[element]) {
-      fPosPedSum[element] += adcpos;
+      fPosPedSum[element]  += adcpos;
       fPosPedSum2[element] += adcpos*adcpos;
       fPosPedCount[element]++;
-      if(fPosPedCount[element] == fMinPeds/5) {
+      if(fPosPedCount[element] == fMinPeds/5)
 	fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
-      }
     }
     if(adcneg <= fNegPedLimit[element]) {
       fNegPedSum[element] += adcneg;
       fNegPedSum2[element] += adcneg*adcneg;
       fNegPedCount[element]++;
-      if(fNegPedCount[element] == fMinPeds/5) {
+      if(fNegPedCount[element] == fMinPeds/5)
 	fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
-      }
     }
     ihit++;
   }
-
   fNPedestalEvents++;
-
   return;
 }
 
 //_____________________________________________________________________________
-void THcAerogel::CalculatePedestals( )
+// Method for calculating pedestals in the 6 GeV era
+void THcAerogel::CalculatePedestals()
 {
   // Use the accumulated pedestal data to calculate pedestals
   // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
@@ -705,13 +925,11 @@ void THcAerogel::CalculatePedestals( )
   for(Int_t i=0; i<fNelem;i++) {
 
     // Positive tubes
-    fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
+    fPosPed[i]    = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
     fPosThresh[i] = fPosPed[i] + 15;
-
     // Negative tubes
-    fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
+    fNegPed[i]    = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
     fNegThresh[i] = fNegPed[i] + 15;
-
     //    cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
 
     // Just a copy for now, but allow the possibility that fXXXPedMean is set
@@ -719,35 +937,34 @@ void THcAerogel::CalculatePedestals( )
     // pedestal events.  (So that pedestals are sensible even if the pedestal events were
     // not acquired.)
     if(fMinPeds > 0) {
-      if(fPosPedCount[i] > fMinPeds) {
+      if(fPosPedCount[i] > fMinPeds)
 	fPosPedMean[i] = fPosPed[i];
-      }
-      if(fNegPedCount[i] > fMinPeds) {
+      if(fNegPedCount[i] > fMinPeds)
 	fNegPedMean[i] = fNegPed[i];
-      }
     }
   }
+}
 
-  //  cout << " " << endl;
-
+//_____________________________________________________________________________
+Int_t THcAerogel::GetIndex(Int_t nRegion, Int_t nValue)
+{
+  return fNRegions * nValue + nRegion;
 }
-void THcAerogel::Print( const Option_t* opt) const {
+
+//_____________________________________________________________________________
+void THcAerogel::Print(const Option_t* opt) const
+{
   THaNonTrackingDetector::Print(opt);
 
   // Print out the pedestals
-
   cout << endl;
   cout << "Aerogel Pedestals" << endl;
   cout << "No.   Neg    Pos" << endl;
-  for(Int_t i=0; i<fNelem; i++){
-    cout << " " << i << "    " << fNegPedMean[i] << "    " << fPosPedMean[i]
-	 << endl;
-  }
+  for(Int_t i=0; i<fNelem; i++)
+    cout << " " << i << "\t" << fNegPedMean[i] << "\t" << fPosPedMean[i] << endl;
   cout << endl;
-
   cout << " fMinPeds = " << fMinPeds << endl;
   cout << endl;
-
 }
 
 ClassImp(THcAerogel)
diff --git a/src/THcAerogel.h b/src/THcAerogel.h
index fb40ebf..1e68df1 100644
--- a/src/THcAerogel.h
+++ b/src/THcAerogel.h
@@ -15,109 +15,154 @@
 class THcAerogel : public THaNonTrackingDetector, public THcHitList {
 
  public:
-  THcAerogel( const char* name, const char* description = "",
-		THaApparatus* a = NULL );
+  THcAerogel(const char* name, const char* description = "", THaApparatus* a = NULL);
   virtual ~THcAerogel();
 
-  virtual void 	     Clear( Option_t* opt="" );
-  virtual Int_t      Decode( const THaEvData& );
-  void               InitArrays();
-  void               DeleteArrays();
-  virtual EStatus    Init( const TDatime& run_time );
-  virtual Int_t      ReadDatabase( const TDatime& date );
-  virtual Int_t      DefineVariables( EMode mode = kDefine );
-  virtual Int_t      CoarseProcess( TClonesArray& tracks );
-  virtual Int_t      FineProcess( TClonesArray& tracks );
-
-  virtual void AccumulatePedestals(TClonesArray* rawhits);
-  virtual void CalculatePedestals();
-
-  virtual Int_t      ApplyCorrections( void );
-
-  virtual void Print(const Option_t* opt) const;
+  virtual void 	  Clear(Option_t* opt="");
+  virtual void    Print(const Option_t* opt) const;
+  virtual void    AccumulatePedestals(TClonesArray* rawhits);
+  virtual void    CalculatePedestals();
+  virtual Int_t   Decode(const THaEvData&);
+  virtual Int_t   ReadDatabase(const TDatime& date);
+  virtual Int_t   DefineVariables(EMode mode = kDefine);
+  virtual Int_t   CoarseProcess(TClonesArray& tracks);
+  virtual Int_t   FineProcess(TClonesArray& tracks);
+  virtual Int_t   ApplyCorrections(void);
+  virtual EStatus Init(const TDatime& run_time);
+
+  void  InitArrays();
+  void  DeleteArrays();
+  Int_t GetIndex(Int_t nRegion, Int_t nValue);
 
   THcAerogel();  // for ROOT I/O
- protected:
-  Int_t fAnalyzePedestals;
 
-  // Parameters
-  Double_t* fPosGain;
-  Double_t* fNegGain;
+ protected:
 
   // Event information
   Int_t fNhits;
 
-  Float_t*   fA_Pos;         // [fNelem] Array of ADC amplitudes
-  Float_t*   fA_Neg;         // [fNelem] Array of ADC amplitudes
-  Float_t*   fA_Pos_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
-  Float_t*   fA_Neg_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
-  Float_t*   fT_Pos;         // [fNelem] Array of TDCs
-  Float_t*   fT_Neg;         // [fNelem] Array of TDCs
-
-  Double_t fPosNpeSum;
-  Double_t fNegNpeSum;
-  Double_t fNpeSum;
-  Int_t fNGoodHits;
-  Int_t fNADCPosHits;
-  Int_t fNADCNegHits;
-  Int_t fNTDCPosHits;
-  Int_t fNTDCNegHits;
-
-  Double_t* fPosNpe;		// [fNelem] # Photoelectrons per positive tube
-  Double_t* fNegNpe;		// [fNelem] # Photoelectrons per negative tube
-
-  // Hits
-  TClonesArray* fPosTDCHits;
-  TClonesArray* fNegTDCHits;
-  TClonesArray* fPosADCHits;
-  TClonesArray* fNegADCHits;
-
-  // Pedestals
-  Int_t fNPedestalEvents;
-  Int_t fMinPeds;
-  Int_t *fPosPedSum;		/* Accumulators for pedestals */
-  Int_t *fPosPedSum2;
-  Int_t *fPosPedLimit;
-  Int_t *fPosPedCount;
-  Int_t *fNegPedSum;
-  Int_t *fNegPedSum2;
-  Int_t *fNegPedLimit;
-  Int_t *fNegPedCount;
-
-  Double_t *fPosPed;
-  Double_t *fPosSig;
-  Double_t *fPosThresh;
-  Double_t *fNegPed;
-  Double_t *fNegSig;
-  Double_t *fNegThresh;
-
-  Double_t *fPosPedMean; 	/* Can be supplied in parameters and then */
-  Double_t *fNegPedMean;	/* be overwritten from ped analysis */
-
-  Int_t fTdcOffset; /* Global TDC offset */
-
+  // 12 GeV variables
+  // Vector/TClonesArray length parameters
+  static const Int_t MaxNumPosAeroPmt = 7;
+  static const Int_t MaxNumNegAeroPmt = 7;
+  static const Int_t MaxNumAdcPulse   = 4;
+  // Tracking variables
+  Int_t     fNRegions;
+  Int_t     fRegionsValueMax;
+  Int_t     fDebugAdc;
+  Double_t  fRedChi2Min;
+  Double_t  fRedChi2Max;
+  Double_t  fBetaMin;
+  Double_t  fBetaMax;
+  Double_t  fENormMin;
+  Double_t  fENormMax;
+  Double_t  fDiffBoxZPos;
+  Double_t  fNpeThresh;
+  Double_t  fAdcTimeWindowMin;
+  Double_t  fAdcTimeWindowMax;
+  Double_t  *fRegionValue;
+  // Counting variables
+  Int_t     fTotNumAdcHits;
+  Int_t     fTotNumGoodAdcHits;
+  Int_t     fTotNumPosAdcHits;
+  Int_t     fTotNumGoodPosAdcHits;
+  Int_t     fTotNumNegAdcHits;
+  Int_t     fTotNumGoodNegAdcHits;
+  Int_t     fTotNumTracksMatched;
+  Int_t     fTotNumTracksFired;
+  // NPE variables
+  Double_t  fPosNpeSum;
+  Double_t  fNegNpeSum;
+  Double_t  fNpeSum;
+  Double_t  *fPosGain;
+  Double_t  *fNegGain;
+  // FADC data objects
   TClonesArray* frPosAdcPedRaw;
   TClonesArray* frPosAdcPulseIntRaw;
   TClonesArray* frPosAdcPulseAmpRaw;
   TClonesArray* frPosAdcPulseTimeRaw;
-
   TClonesArray* frPosAdcPed;
   TClonesArray* frPosAdcPulseInt;
   TClonesArray* frPosAdcPulseAmp;
-
   TClonesArray* frNegAdcPedRaw;
   TClonesArray* frNegAdcPulseIntRaw;
   TClonesArray* frNegAdcPulseAmpRaw;
   TClonesArray* frNegAdcPulseTimeRaw;
-
   TClonesArray* frNegAdcPed;
   TClonesArray* frNegAdcPulseInt;
   TClonesArray* frNegAdcPulseAmp;
+  TClonesArray* fPosAdcErrorFlag;
+  TClonesArray* fNegAdcErrorFlag;
+  // Individual PMT data objects
+  vector<Int_t>    fNumPosAdcHits;
+  vector<Int_t>    fNumNegAdcHits;
+  vector<Int_t>    fNumGoodPosAdcHits;
+  vector<Int_t>    fNumGoodNegAdcHits;
+  vector<Int_t>    fNumTracksMatched;
+  vector<Int_t>    fNumTracksFired;
+  vector<Double_t> fPosNpe;
+  vector<Double_t> fNegNpe;
+  vector<Double_t> fGoodPosAdcPed;
+  vector<Double_t> fGoodPosAdcPulseInt;
+  vector<Double_t> fGoodPosAdcPulseIntRaw;
+  vector<Double_t> fGoodPosAdcPulseAmp;
+  vector<Double_t> fGoodPosAdcPulseTime;
+  vector<Double_t> fGoodNegAdcPed;
+  vector<Double_t> fGoodNegAdcPulseInt;
+  vector<Double_t> fGoodNegAdcPulseIntRaw;
+  vector<Double_t> fGoodNegAdcPulseAmp;
+  vector<Double_t> fGoodNegAdcPulseTime;
+
+  // 6 GeV era variables
+  Int_t     fAnalyzePedestals;
+  Int_t     fSixGevData;
+  Int_t     fNGoodHits;
+  Int_t     fNADCPosHits;
+  Int_t     fNADCNegHits;
+  Int_t     fNTDCPosHits;
+  Int_t     fNTDCNegHits;
+  Int_t     fTdcOffset; /* Global TDC offset */
+  Int_t     fNPedestalEvents;
+  Int_t     fMinPeds;
+  Double_t  fPosNpeSumSixGev;
+  Double_t  fNegNpeSumSixGev;
+  Double_t  fNpeSumSixGev;
+  Int_t    *fPosPedSum;		/* Accumulators for pedestals */
+  Int_t    *fPosPedSum2;
+  Int_t    *fPosPedLimit;
+  Int_t    *fPosPedCount;
+  Int_t    *fNegPedSum;
+  Int_t    *fNegPedSum2;
+  Int_t    *fNegPedLimit;
+  Int_t    *fNegPedCount;
+  Float_t  *fA_Pos;          // [fNelem] Array of ADC amplitudes
+  Float_t  *fA_Neg;          // [fNelem] Array of ADC amplitudes
+  Float_t  *fA_Pos_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
+  Float_t  *fA_Neg_p;	     // [fNelem] Array of ped-subtracted ADC amplitudes
+  Float_t  *fT_Pos;          // [fNelem] Array of TDCs
+  Float_t  *fT_Neg;          // [fNelem] Array of TDCs
+  Double_t *fPosPed;
+  Double_t *fPosSig;
+  Double_t *fPosThresh;
+  Double_t *fNegPed;
+  Double_t *fNegSig;
+  Double_t *fNegThresh;
+  Double_t *fPosPedMean; 	/* Can be supplied in parameters and then */
+  Double_t *fNegPedMean;	/* be overwritten from ped analysis */
+
+  TClonesArray *fPosTDCHits;
+  TClonesArray *fNegTDCHits;
+  TClonesArray *fPosADCHits;
+  TClonesArray *fNegADCHits;
+
+  vector<Double_t> fPosNpeSixGev;
+  vector<Double_t> fNegNpeSixGev;
 
   void Setup(const char* name, const char* description);
   virtual void  InitializePedestals( );
 
   ClassDef(THcAerogel,0)   // Generic aerogel class
-};
+}
+;
 
 #endif
diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx
index ceb1aef..cea6b3d 100644
--- a/src/THcCherenkov.cxx
+++ b/src/THcCherenkov.cxx
@@ -1,9 +1,7 @@
 /** \class THcCherenkov
     \ingroup Detectors
 
-Class for an Cherenkov detector consisting of two PMT's
-
-\author Zafar Ahmed
+    Class for Cherenkov detectors
 
 */
 
@@ -23,7 +21,6 @@ Class for an Cherenkov detector consisting of two PMT's
 #include "THaTrack.h"
 #include "TClonesArray.h"
 #include "TMath.h"
-
 #include "THaTrackProj.h"
 
 #include <algorithm>
@@ -32,14 +29,13 @@ Class for an Cherenkov detector consisting of two PMT's
 #include <cstdlib>
 #include <iostream>
 #include <string>
+#include <iomanip>
 
 using namespace std;
 
 using std::cout;
 using std::cin;
 using std::endl;
-
-#include <iomanip>
 using std::setw;
 using std::setprecision;
 
@@ -49,20 +45,27 @@ THcCherenkov::THcCherenkov( const char* name, const char* description,
   THaNonTrackingDetector(name,description,apparatus)
 {
   // Normal constructor with name and description
-  fADCHits = new TClonesArray("THcSignalHit",16);
-
-  frAdcPedRaw = new TClonesArray("THcSignalHit", 16);
-  frAdcPulseIntRaw = new TClonesArray("THcSignalHit", 16);
-  frAdcPulseAmpRaw = new TClonesArray("THcSignalHit", 16);
-  frAdcPulseTimeRaw = new TClonesArray("THcSignalHit", 16);
-
-  frAdcPed = new TClonesArray("THcSignalHit", 16);
-  frAdcPulseInt = new TClonesArray("THcSignalHit", 16);
-  frAdcPulseAmp = new TClonesArray("THcSignalHit", 16);
+  frAdcPedRaw       = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPulseIntRaw  = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPulseAmpRaw  = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPulseTimeRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPed          = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPulseInt     = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  frAdcPulseAmp     = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+  fAdcErrorFlag     = new TClonesArray("THcSignalHit", MaxNumCerPmt*MaxNumAdcPulse);
+
+  fNumAdcHits         = vector<Int_t>    (MaxNumCerPmt, 0.0);
+  fNumGoodAdcHits     = vector<Int_t>    (MaxNumCerPmt, 0.0);
+  fNumTracksMatched   = vector<Int_t>    (MaxNumCerPmt, 0.0);
+  fNumTracksFired     = vector<Int_t>    (MaxNumCerPmt, 0.0);
+  fNpe                = vector<Double_t> (MaxNumCerPmt, 0.0);
+  fGoodAdcPed         = vector<Double_t> (MaxNumCerPmt, 0.0);
+  fGoodAdcPulseInt    = vector<Double_t> (MaxNumCerPmt, 0.0);
+  fGoodAdcPulseIntRaw = vector<Double_t> (MaxNumCerPmt, 0.0);
+  fGoodAdcPulseAmp    = vector<Double_t> (MaxNumCerPmt, 0.0);
+  fGoodAdcPulseTime   = vector<Double_t> (MaxNumCerPmt, 0.0);
 
   InitArrays();
-
-
 }
 
 //_____________________________________________________________________________
@@ -70,81 +73,69 @@ THcCherenkov::THcCherenkov( ) :
   THaNonTrackingDetector()
 {
   // Constructor
-  fADCHits = NULL;
-
-  frAdcPedRaw = NULL;
-  frAdcPulseIntRaw = NULL;
-  frAdcPulseAmpRaw = NULL;
+  frAdcPedRaw       = NULL;
+  frAdcPulseIntRaw  = NULL;
+  frAdcPulseAmpRaw  = NULL;
   frAdcPulseTimeRaw = NULL;
-
-  frAdcPed = NULL;
-  frAdcPulseInt = NULL;
-  frAdcPulseAmp = NULL;
+  frAdcPed          = NULL;
+  frAdcPulseInt     = NULL;
+  frAdcPulseAmp     = NULL;
+  fAdcErrorFlag     = NULL;
 
   InitArrays();
 }
 
+//_____________________________________________________________________________
+THcCherenkov::~THcCherenkov()
+{
+  // Destructor
+  delete frAdcPedRaw;       frAdcPedRaw       = NULL;
+  delete frAdcPulseIntRaw;  frAdcPulseIntRaw  = NULL;
+  delete frAdcPulseAmpRaw;  frAdcPulseAmpRaw  = NULL;
+  delete frAdcPulseTimeRaw; frAdcPulseTimeRaw = NULL;
+  delete frAdcPed;          frAdcPed          = NULL;
+  delete frAdcPulseInt;     frAdcPulseInt     = NULL;
+  delete frAdcPulseAmp;     frAdcPulseAmp     = NULL;
+  delete fAdcErrorFlag;     fAdcErrorFlag     = NULL;
+
+  DeleteArrays();
+}
+
 //_____________________________________________________________________________
 void THcCherenkov::InitArrays()
 {
-  fGain = NULL;
-  fCerWidth = NULL;
-  fNPMT = NULL;
-  fADC = NULL;
-  fADC_P = NULL;
-  fNPE = NULL;
-  fPedSum = NULL;
-  fPedSum2 = NULL;
+  fGain     = NULL;
+  fPedSum   = NULL;
+  fPedSum2  = NULL;
   fPedLimit = NULL;
-  fPedMean = NULL;
+  fPedMean  = NULL;
   fPedCount = NULL;
-  fPed = NULL;
-  fThresh = NULL;
+  fPed      = NULL;
+  fThresh   = NULL;
 }
 //_____________________________________________________________________________
 void THcCherenkov::DeleteArrays()
 {
   delete [] fGain; fGain = NULL;
-  delete [] fCerWidth; fCerWidth = NULL;
-  delete [] fNPMT; fNPMT = NULL;
-  delete [] fADC; fADC = NULL;
-  delete [] fADC; fADC_P = NULL;
-  delete [] fNPE; fNPE = NULL;
-  delete [] fPedSum; fPedSum = NULL;
-  delete [] fPedSum2; fPedSum2 = NULL;
+
+  // 6 Gev variables
+  delete [] fPedSum;   fPedSum   = NULL;
+  delete [] fPedSum2;  fPedSum2  = NULL;
   delete [] fPedLimit; fPedLimit = NULL;
-  delete [] fPedMean; fPedMean = NULL;
+  delete [] fPedMean;  fPedMean  = NULL;
   delete [] fPedCount; fPedCount = NULL;
-  delete [] fPed; fPed = NULL;
-  delete [] fThresh; fThresh = NULL;
-}
-//_____________________________________________________________________________
-THcCherenkov::~THcCherenkov()
-{
-  // Destructor
-  delete fADCHits; fADCHits = NULL;
-
-  delete frAdcPedRaw; frAdcPedRaw = NULL;
-  delete frAdcPulseIntRaw; frAdcPulseIntRaw = NULL;
-  delete frAdcPulseAmpRaw; frAdcPulseAmpRaw = NULL;
-  delete frAdcPulseTimeRaw; frAdcPulseTimeRaw = NULL;
-
-  delete frAdcPed; frAdcPed = NULL;
-  delete frAdcPulseInt; frAdcPulseInt = NULL;
-  delete frAdcPulseAmp; frAdcPulseAmp = NULL;
-
-  DeleteArrays();
-
+  delete [] fPed;      fPed      = NULL;
+  delete [] fThresh;   fThresh   = NULL;
 }
 
 //_____________________________________________________________________________
 THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
 {
-  cout << "THcCherenkov::Init " << GetName() << endl;
+  cout << "THcCherenkov::Init for: " << GetName() << endl;
 
   string EngineDID = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
   std::transform(EngineDID.begin(), EngineDID.end(), EngineDID.begin(), ::toupper);
-  if( gHcDetectorMap->FillMap(fDetMap, EngineDID.c_str()) < 0 ) {
+  if(gHcDetectorMap->FillMap(fDetMap, EngineDID.c_str()) < 0) {
     static const char* const here = "Init()";
     Error(Here(here), "Error filling detectormap for %s.", EngineDID.c_str());
     return kInitError;
@@ -155,7 +146,7 @@ THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
   InitHitList(fDetMap, "THcCherenkovHit", fDetMap->GetTotNumChan()+1);
 
   EStatus status;
-  if( (status = THaNonTrackingDetector::Init( date )) )
+  if((status = THaNonTrackingDetector::Init( date )))
     return fStatus=status;
 
   return fStatus = kOK;
@@ -167,70 +158,67 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date )
   // This function is called by THaDetectorBase::Init() once at the beginning
   // of the analysis.
 
-  cout << "THcCherenkov::ReadDatabase " << GetName() << endl; // Ahmed
+  cout << "THcCherenkov::ReadDatabase for: " << GetName() << endl; // Ahmed
 
   string prefix = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
   std::transform(prefix.begin(), prefix.end(), prefix.begin(), ::tolower);
 
+  fNRegions = 4;  // Defualt if not set in paramter file
+
   DBRequest list_1[] = {
     {"_tot_pmts", &fNelem, kInt},
     {0}
   };
+
   gHcParms->LoadParmValues(list_1, prefix.c_str());
 
-  //    fNelem = 2;      // Default if not defined
-  fCerNRegions = 3;
+  Bool_t optional = true;
 
-  fNPMT = new Int_t[fNelem];
-  fADC_hit = new Int_t[fNelem];
-  fADC = new Double_t[fNelem];
-  fADC_P = new Double_t[fNelem];
-  fNPE = new Double_t[fNelem];
+  cout << "Number of " << GetApparatus()->GetName() << "."
+       << GetName() << " PMTs defined = " << fNelem << endl;
 
-  fCerWidth = new Double_t[fNelem];
-  fGain = new Double_t[fNelem];
+  // 6 GeV pedestal paramters
   fPedLimit = new Int_t[fNelem];
-  fPedMean = new Double_t[fNelem];
+  fGain     = new Double_t[fNelem];
+  fPedMean  = new Double_t[fNelem];
 
- 
-
-  fCerTrackCounter = new Int_t [fCerNRegions];
-  fCerFiredCounter = new Int_t [fCerNRegions];
-  for ( Int_t ireg = 0; ireg < fCerNRegions; ireg++ ) {
-    fCerTrackCounter[ireg] = 0;
-    fCerFiredCounter[ireg] = 0;
-  }
-
-  fCerRegionsValueMax = fCerNRegions * 8; // This value 8 should also be in paramter file
-  fCerRegionValue = new Double_t [fCerRegionsValueMax];
+  // Region parameters
+  fRegionsValueMax = fNRegions * 8;
+  fRegionValue     = new Double_t[fRegionsValueMax];
 
   DBRequest list[]={
-    {"_adc_to_npe", fGain,     kDouble, (UInt_t) fNelem},
-    {"_ped_limit",  fPedLimit, kInt,    (UInt_t) fNelem},
-    {"_width",      fCerWidth, kDouble, (UInt_t) fNelem},
-    {"_chi2max",     &fCerChi2Max,        kDouble},
-    {"_beta_min",    &fCerBetaMin,        kDouble},
-    {"_beta_max",    &fCerBetaMax,        kDouble},
-    {"_et_min",      &fCerETMin,          kDouble},
-    {"_et_max",      &fCerETMax,          kDouble},
-    {"_mirror_zpos", &fCerMirrorZPos,     kDouble},
-    {"_region",      &fCerRegionValue[0], kDouble, (UInt_t) fCerRegionsValueMax},
-    {"_threshold",   &fCerThresh,         kDouble},
-    //    {"cer_regions",     &fCerNRegions,       kInt},
+    {"_ped_limit",        fPedLimit,          kInt,     (UInt_t) fNelem, optional},
+    {"_adc_to_npe",       fGain,              kDouble,  (UInt_t) fNelem},
+    {"_red_chi2_min",     &fRedChi2Min,       kDouble},
+    {"_red_chi2_max",     &fRedChi2Max,       kDouble},
+    {"_beta_min",         &fBetaMin,          kDouble},
+    {"_beta_max",         &fBetaMax,          kDouble},
+    {"_enorm_min",        &fENormMin,         kDouble},
+    {"_enorm_max",        &fENormMax,         kDouble},
+    {"_mirror_zpos",      &fMirrorZPos,       kDouble},
+    {"_npe_thresh",       &fNpeThresh,        kDouble},
+    {"_debug_adc",        &fDebugAdc,         kInt, 0, 1},
+    {"_adcTimeWindowMin", &fAdcTimeWindowMin, kDouble},
+    {"_adcTimeWindowMax", &fAdcTimeWindowMax, kDouble},
+    {"_num_regions",      &fNRegions,         kInt},
+    {"_region",           &fRegionValue[0],   kDouble,  (UInt_t) fRegionsValueMax},
     {0}
   };
 
-  gHcParms->LoadParmValues((DBRequest*)&list,prefix.c_str());
+  fDebugAdc = 0; // Set ADC debug parameter to false unless set in parameter file
 
-  fIsInit = true;
+  gHcParms->LoadParmValues((DBRequest*)&list, prefix.c_str());
 
+  if (fDebugAdc) cout << "Cherenkov ADC Debug Flag Set To TRUE" << endl;
 
-  for (Int_t i1 = 0; i1 < fCerNRegions; i1++ ) {
-    cout << "Region " << i1 << endl;
-    for (Int_t i2 = 0; i2 < 8; i2++ ) {
-      cout << fCerRegionValue[GetCerIndex( i1, i2 )] << " ";
-    }
-    cout <<endl;
+  fIsInit = true;
+
+  cout << "Track Matching Parameters for: " << GetName() << endl;
+  for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
+    cout << "Region = " << iregion + 1 << endl;
+    for (Int_t ivalue = 0; ivalue < 8; ivalue++)
+      cout << fRegionValue[GetIndex(iregion, ivalue)] << "  ";
+    cout << endl;
   }
 
   // Create arrays to hold pedestal results
@@ -243,64 +231,64 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date )
 Int_t THcCherenkov::DefineVariables( EMode mode )
 {
   // Initialize global variables for histogramming and tree
-
-  cout << "THcCherenkov::DefineVariables called " << GetName() << endl;
+  cout << "THcCherenkov::DefineVariables called for: " << GetName() << endl;
 
   if( mode == kDefine && fIsSetup ) return kOK;
   fIsSetup = ( mode == kDefine );
 
   // Register variables in global list
+  vector<RVarDef> vars;
+
+  vars.push_back({"adcCounter",   "ADC counter numbers",            "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"});
+  vars.push_back({"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"});
+
+  vars.push_back({"numGoodAdcHits",    "Number of Good ADC Hits Per PMT", "fNumGoodAdcHits"});    // Cherenkov occupancy
+  vars.push_back({"totNumGoodAdcHits", "Total Number of Good ADC Hits",   "fTotNumGoodAdcHits"}); // Cherenkov multiplicity
+
+  vars.push_back({"numTracksMatched",    "Number of Tracks Matched Per Region",       "fNumTracksMatched"});
+  vars.push_back({"numTracksFired",      "Number of Tracks that Fired Per Region",    "fNumTracksFired"});
+  vars.push_back({"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"});
+  vars.push_back({"totNumTracksFired",   "Total Number of Tracks that Fired",         "fTotNumTracksFired"});
+
+  vars.push_back({"npe",          "Number of PEs",                  "fNpe"});
+  vars.push_back({"npeSum",       "Total Number of PEs",            "fNpeSum"});
+
+  vars.push_back({"goodAdcPed",          "Good ADC pedestals",           "fGoodAdcPed"});
+  vars.push_back({"goodAdcPulseInt",     "Good ADC pulse integrals",     "fGoodAdcPulseInt"});
+  vars.push_back({"goodAdcPulseIntRaw",  "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"});
+  vars.push_back({"goodAdcPulseAmp",     "Good ADC pulse amplitudes",    "fGoodAdcPulseAmp"});
+  vars.push_back({"goodAdcPulseTime",    "Good ADC pulse times",         "fGoodAdcPulseTime"});
+
+  if (fDebugAdc) {
+    vars.push_back({"numAdcHits",      "Number of ADC Hits Per PMT", "fNumAdcHits"});        // Cherenkov occupancy
+    vars.push_back({"totNumAdcHits",   "Total Number of ADC Hits",   "fTotNumAdcHits"});     // Cherenkov multiplicity
+    vars.push_back({"adcPedRaw",       "Raw ADC pedestals",          "frAdcPedRaw.THcSignalHit.GetData()"});
+    vars.push_back({"adcPulseIntRaw",  "Raw ADC pulse integrals",    "frAdcPulseIntRaw.THcSignalHit.GetData()"});
+    vars.push_back({"adcPulseAmpRaw",  "Raw ADC pulse amplitudes",   "frAdcPulseAmpRaw.THcSignalHit.GetData()"});
+    vars.push_back({"adcPulseTimeRaw", "Raw ADC pulse times",        "frAdcPulseTimeRaw.THcSignalHit.GetData()"});
+    vars.push_back({"adcPed",          "ADC pedestals",              "frAdcPed.THcSignalHit.GetData()"});
+    vars.push_back({"adcPulseInt",     "ADC pulse integrals",        "frAdcPulseInt.THcSignalHit.GetData()"});
+    vars.push_back({"adcPulseAmp",     "ADC pulse amplitudes",       "frAdcPulseAmp.THcSignalHit.GetData()"});
+  }
 
-  // Do we need to put the number of pos/neg TDC/ADC hits into the variables?
-  // No.  They show up in tree as Ndata.H.aero.postdchits for example
-
-  RVarDef vars[] = {
-    {"phototubes",      "Nuber of Cherenkov photo tubes",        "fNPMT"},
-    {"adc",             "Raw ADC values",                        "fADC"},
-    {"adc_hit",         "ADC hit flag =1 means hit",             "fADC_hit"},
-    {"adc_p",           "Pedestal Subtracted ADC values",        "fADC_P"},
-    {"npe",             "Number of Photo electrons",             "fNPE"},
-    {"npesum",          "Sum of Number of Photo electrons",      "fNPEsum"},
-    {"ncherhit",        "Number of Hits(Cherenkov)",             "fNCherHit"},
-    {"certrackcounter", "Tracks inside Cherenkov region",        "fCerTrackCounter"},
-    {"cerfiredcounter", "Tracks with engough Cherenkov NPEs ",   "fCerFiredCounter"},
-
-    {"adcCounter",      "List of ADC counter numbers.",      "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
-
-    {"adcPedRaw",       "List of raw ADC pedestals",         "frAdcPedRaw.THcSignalHit.GetData()"},
-    {"adcPulseIntRaw",  "List of raw ADC pulse integrals.",  "frAdcPulseIntRaw.THcSignalHit.GetData()"},
-    {"adcPulseAmpRaw",  "List of raw ADC pulse amplitudes.", "frAdcPulseAmpRaw.THcSignalHit.GetData()"},
-    {"adcPulseTimeRaw", "List of raw ADC pulse times.",      "frAdcPulseTimeRaw.THcSignalHit.GetData()"},
+  RVarDef end {};
+  vars.push_back(end);
 
-    {"adcPed",          "List of ADC pedestals",             "frAdcPed.THcSignalHit.GetData()"},
-    {"adcPulseInt",     "List of ADC pulse integrals.",      "frAdcPulseInt.THcSignalHit.GetData()"},
-    {"adcPulseAmp",     "List of ADC pulse amplitudes.",     "frAdcPulseAmp.THcSignalHit.GetData()"},
+  return DefineVarsFromList(vars.data(), mode);
 
-    { 0 }
-  };
-
-  return DefineVarsFromList( vars, mode );
 }
 //_____________________________________________________________________________
 inline
 void THcCherenkov::Clear(Option_t* opt)
 {
   // Clear the hit lists
-  fADCHits->Clear();
-
-  // Clear Cherenkov variables  from h_trans_cer.f
+  fNhits               = 0;
+  fTotNumAdcHits       = 0;
+  fTotNumGoodAdcHits   = 0;
+  fTotNumTracksMatched = 0;
+  fTotNumTracksFired   = 0;
 
-  fNhits = 0;	     
-  fNPEsum = 0.0;
-  fNCherHit = 0;
-
-  for(Int_t itube = 0;itube < fNelem;itube++) {
-    fNPMT[itube] = 0;
-    fADC[itube] = 0;
-    fADC_P[itube] = 0;
-    fNPE[itube] = 0;
-    fADC_hit[itube] = 0;
- }
+  fNpeSum = 0.0;
 
   frAdcPedRaw->Clear();
   frAdcPulseIntRaw->Clear();
@@ -310,6 +298,24 @@ void THcCherenkov::Clear(Option_t* opt)
   frAdcPed->Clear();
   frAdcPulseInt->Clear();
   frAdcPulseAmp->Clear();
+  fAdcErrorFlag->Clear();
+
+  for (UInt_t ielem = 0; ielem < fNumAdcHits.size(); ielem++)
+    fNumAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumGoodAdcHits.size(); ielem++)
+    fNumGoodAdcHits.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumTracksMatched.size(); ielem++)
+    fNumTracksMatched.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fNumTracksFired.size(); ielem++)
+    fNumTracksFired.at(ielem) = 0;
+  for (UInt_t ielem = 0; ielem < fGoodAdcPed.size(); ielem++) {
+    fGoodAdcPed.at(ielem)         = 0.0;
+    fGoodAdcPulseInt.at(ielem)    = 0.0;
+    fGoodAdcPulseIntRaw.at(ielem) = 0.0;
+    fGoodAdcPulseAmp.at(ielem)    = 0.0;
+    fGoodAdcPulseTime.at(ielem)   = 0.0;
+    fNpe.at(ielem)                = 0.0;
+  }
 
 }
 
@@ -330,40 +336,35 @@ Int_t THcCherenkov::Decode( const THaEvData& evdata )
     fAnalyzePedestals = 0;	// Don't analyze pedestals next event
   }
 
-
-  Int_t ihit = 0;
-  Int_t nADCHits=0;
-
+  Int_t  ihit      = 0;
   UInt_t nrAdcHits = 0;
 
   while(ihit < fNhits) {
-    THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit);
 
-    Int_t padnum = hit->fCounter;
+    THcCherenkovHit* hit       = (THcCherenkovHit*) fRawHitList->At(ihit);
+    Int_t            npmt      = hit->fCounter;
+    THcRawAdcHit&    rawAdcHit = hit->GetRawAdcHitPos();
 
-    THcRawAdcHit& rawAdcHit = hit->GetRawAdcHitPos();
-    for (UInt_t thit=0; thit<rawAdcHit.GetNPulses(); ++thit) {
-      ((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPedRaw());
-      ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPed());
+    for (UInt_t thit = 0; thit < rawAdcHit.GetNPulses(); thit++) {
 
-      ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseIntRaw(thit));
-      ((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseInt(thit));
+      ((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPedRaw());
+      ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPed());
 
-      ((THcSignalHit*) frAdcPulseAmpRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseAmpRaw(thit));
-      ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseAmp(thit));
+      ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPulseIntRaw(thit));
+      ((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPulseInt(thit));
 
-      ((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(nrAdcHits))->Set(padnum, rawAdcHit.GetPulseTimeRaw(thit));
-      fADC_hit[padnum-1]=1; // 
-      //     cout << dec << "thit = " << thit << " " << padnum << " " << rawAdcHit.GetPulseInt(thit)<< " " << rawAdcHit.GetPulseIntRaw(thit) << " " << rawAdcHit.GetPedRaw() << " " << rawAdcHit.GetPedRaw()*28./4.<< endl;
-      ++nrAdcHits;
-    }
+      ((THcSignalHit*) frAdcPulseAmpRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPulseAmpRaw(thit));
+      ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPulseAmp(thit));
 
-    // ADC hit
-    if(hit->GetRawAdcHitPos().GetPulseIntRaw() >  0) {
-      THcSignalHit *sighit = (THcSignalHit*) fADCHits->ConstructedAt(nADCHits++);
-      sighit->Set(hit->fCounter, hit->GetRawAdcHitPos().GetPulseIntRaw());
-    }
+      ((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPulseTimeRaw(thit));
 
+      if (rawAdcHit.GetPulseAmpRaw(thit) > 0)  ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 0);
+      if (rawAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 1);
+
+      ++nrAdcHits;
+      fTotNumAdcHits++;
+      fNumAdcHits.at(npmt-1) = npmt;
+    }
     ihit++;
   }
   return ihit;
@@ -376,24 +377,36 @@ Int_t THcCherenkov::ApplyCorrections( void )
 }
 
 //_____________________________________________________________________________
-Int_t THcCherenkov::CoarseProcess( TClonesArray&  ) //tracks
+Int_t THcCherenkov::CoarseProcess( TClonesArray&  )
 {
-  for(Int_t ihit=0; ihit < fNhits; ihit++) {
-    THcCherenkovHit* hit = (THcCherenkovHit *) fRawHitList->At(ihit); // nhit = 1, hcer_tot_hits
-
-    // Pedestal subtraction and gain adjustment
-
-    Int_t npmt = hit->fCounter - 1;                             // tube = hcer_tube_num(nhit)
-    fNPMT[npmt] = hit->fCounter;
-    fADC[npmt] = hit->GetRawAdcHitPos().GetPulseIntRaw();
-    fADC_P[npmt] = hit->GetRawAdcHitPos().GetPulseInt();
 
-    fNPE[npmt] = fGain[npmt]*fADC_P[npmt];
-    fNCherHit ++;
-    fNPEsum += fNPE[npmt];
+  // Loop over the elements in the TClonesArray
+  for(Int_t ielem = 0; ielem < frAdcPulseInt->GetEntries(); ielem++) {
+
+    Int_t    npmt         = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
+    Double_t pulsePed     = ((THcSignalHit*) frAdcPed->ConstructedAt(ielem))->GetData();
+    Double_t pulseInt     = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetData();
+    Double_t pulseIntRaw  = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
+    Double_t pulseAmp     = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
+    Double_t pulseTime    = ((THcSignalHit*) frAdcPulseTimeRaw->ConstructedAt(ielem))->GetData();
+    Bool_t   errorFlag    = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
+    Bool_t   pulseTimeCut = pulseTime > fAdcTimeWindowMin && pulseTime < fAdcTimeWindowMax;
+
+    // By default, the last hit within the timing cut will be considered "good"
+    if (!errorFlag && pulseTimeCut) {
+      fGoodAdcPed.at(npmt)         = pulsePed;
+      fGoodAdcPulseInt.at(npmt)    = pulseInt;
+      fGoodAdcPulseIntRaw.at(npmt) = pulseIntRaw;
+      fGoodAdcPulseAmp.at(npmt)    = pulseAmp;
+      fGoodAdcPulseTime.at(npmt)   = pulseTime;
+
+      fNpe.at(npmt) = fGain[npmt]*fGoodAdcPulseInt.at(npmt);
+      fNpeSum += fNpe.at(npmt);
+
+      fTotNumGoodAdcHits++;
+      fNumGoodAdcHits.at(npmt) = npmt + 1;
+    }
   }
-
-
   return 0;
 }
 
@@ -401,69 +414,82 @@ Int_t THcCherenkov::CoarseProcess( TClonesArray&  ) //tracks
 Int_t THcCherenkov::FineProcess( TClonesArray& tracks )
 {
 
-  if ( tracks.GetLast() > -1 ) {
-
-    THaTrack* theTrack = dynamic_cast<THaTrack*>( tracks.At(0) );
-    if (!theTrack) return -1;
-
-    if ( ( ( tracks.GetLast() + 1 ) == 1 ) &&
-	 ( theTrack->GetChi2()/theTrack->GetNDoF() > 0. ) &&
-	 ( theTrack->GetChi2()/theTrack->GetNDoF() <  fCerChi2Max ) &&
-	 ( theTrack->GetBeta() > fCerBetaMin ) &&
-	 ( theTrack->GetBeta() < fCerBetaMax ) &&
-	 ( ( theTrack->GetEnergy() / theTrack->GetP() ) > fCerETMin ) &&
-	 ( ( theTrack->GetEnergy() / theTrack->GetP() ) < fCerETMax )
-	 ) {
-
-      Double_t cerX = theTrack->GetX() + theTrack->GetTheta() * fCerMirrorZPos;
-      Double_t cerY = theTrack->GetY() + theTrack->GetPhi()   * fCerMirrorZPos;
-
-      for ( Int_t ir = 0; ir < fCerNRegions; ir++ ) {
-
-	//	*     hit must be inside the region in order to continue.
-
-	if ( ( TMath::Abs( fCerRegionValue[GetCerIndex( ir, 0 )] - cerX ) <
-	       fCerRegionValue[GetCerIndex( ir, 4 )] ) &&
-	     ( TMath::Abs( fCerRegionValue[GetCerIndex( ir, 1 )] - cerY ) <
-	       fCerRegionValue[GetCerIndex( ir, 5 )] ) &&
-	     ( TMath::Abs( fCerRegionValue[GetCerIndex( ir, 2 )] - theTrack->GetTheta() ) <
-	       fCerRegionValue[GetCerIndex( ir, 6 )] ) &&
-	     ( TMath::Abs( fCerRegionValue[GetCerIndex( ir, 3 )] - theTrack->GetPhi() ) <
-	       fCerRegionValue[GetCerIndex( ir, 7 )] )
-	     ) {
-
-	  // *     increment the 'should have fired' counters
-	  fCerTrackCounter[ir] ++;
-
-	  // *     increment the 'did fire' counters
-	  if ( fNPEsum > fCerThresh ) {
-	    fCerFiredCounter[ir] ++;
-	  }
-	}
-      } // loop over regions
-    }
-  }
+  Int_t nTracks = tracks.GetLast() + 1;
+
+  for (Int_t itrack = 0; itrack < nTracks; itrack++) {
+
+    THaTrack* track = dynamic_cast<THaTrack*> (tracks[itrack]);
+    if (track->GetIndex() != 0) continue;  // Select the best track
+
+    Double_t trackChi2    = track->GetChi2();
+    Int_t    trackNDoF    = track->GetNDoF();
+    Double_t trackRedChi2 = trackChi2/trackNDoF;
+    Double_t trackBeta    = track->GetBeta();
+    Double_t trackEnergy  = track->GetEnergy();
+    Double_t trackMom     = track->GetP();
+    Double_t trackENorm   = trackEnergy/trackMom;
+    Double_t trackXfp     = track->GetX();
+    Double_t trackYfp     = track->GetY();
+    Double_t trackTheta   = track->GetTheta();
+    Double_t trackPhi     = track->GetPhi();
+
+    Bool_t trackRedChi2Cut = trackRedChi2 > fRedChi2Min && trackRedChi2 < fRedChi2Max;
+    Bool_t trackBetaCut    = trackBeta    > fBetaMin    && trackBeta    < fBetaMax;
+    Bool_t trackENormCut   = trackENorm   > fENormMin   && trackENorm   < fENormMax;
+
+    if (trackRedChi2Cut && trackBetaCut && trackENormCut) {
+
+      // Project the track to the Cherenkov mirror planes
+      Double_t xAtCher = trackXfp + trackTheta * fMirrorZPos;
+      Double_t yAtCher = trackYfp + trackPhi   * fMirrorZPos;
+
+      // cout << "Cherenkov Detector: " << GetName() << " has fNRegions = " << fNRegions << endl;
+      // cout << "nTracks = " << nTracks << "\t" << "trackChi2 = " << trackChi2
+      // 	   << "\t" << "trackNDof = " << trackNDoF << "\t" << "trackRedChi2 = " << trackRedChi2 << endl;
+      // cout << "trackBeta = " << trackBeta << "\t" << "trackEnergy = " << trackEnergy << "\t"
+      // 	   << "trackMom = " << trackMom << "\t" << "trackENorm = " << trackENorm << endl;
+      // cout << "trackXfp = " << trackXfp << "\t" << "trackYfp = " << trackYfp << "\t"
+      // 	   << "trackTheta = " << trackTheta << "\t" << "trackPhi = " << trackPhi << endl;
+      // cout << "fMirrorZPos = " << fMirrorZPos << "\t" << "xAtCher = " << xAtCher << "\t" << "yAtCher = " << yAtCher << endl;
+      // cout << "=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:" << endl;
+
+      for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
+
+	if ((TMath::Abs(fRegionValue[GetIndex(iregion, 0)] - xAtCher)    < fRegionValue[GetIndex(iregion, 4)]) &&
+	    (TMath::Abs(fRegionValue[GetIndex(iregion, 1)] - yAtCher)    < fRegionValue[GetIndex(iregion, 5)]) &&
+	    (TMath::Abs(fRegionValue[GetIndex(iregion, 2)] - trackTheta) < fRegionValue[GetIndex(iregion, 6)]) &&
+	    (TMath::Abs(fRegionValue[GetIndex(iregion, 3)] - trackPhi)   < fRegionValue[GetIndex(iregion, 7)])) {
+
+	  fTotNumTracksMatched++;
+	  fNumTracksMatched.at(iregion) = iregion + 1;
+
+	  if (fNpeSum > fNpeThresh) {
+	    fTotNumTracksFired++;
+	    fNumTracksFired.at(iregion) = iregion + 1;
+	  }  // NPE threshold cut
+	}  // Regional cuts
+      }  // Loop over regions
+    }  // Tracking cuts
+  }  // Track loop
 
   return 0;
 }
 
 //_____________________________________________________________________________
-void THcCherenkov::InitializePedestals( )
+void THcCherenkov::InitializePedestals()
 {
   fNPedestalEvents = 0;
-  fMinPeds = 500; 		// In engine, this is set in parameter file
-  fPedSum = new Int_t [fNelem];
-  fPedSum2 = new Int_t [fNelem];
-  fPedCount = new Int_t [fNelem];
-
-  fPed = new Double_t [fNelem];
-  fThresh = new Double_t [fNelem];
-  for(Int_t i=0;i<fNelem;i++) {
-    fPedSum[i] = 0;
-    fPedSum2[i] = 0;
+  fMinPeds         = 500; 		// In engine, this is set in parameter file
+  fPedSum          = new Int_t [fNelem];
+  fPedSum2         = new Int_t [fNelem];
+  fPedCount        = new Int_t [fNelem];
+  fPed             = new Double_t [fNelem];
+  fThresh          = new Double_t [fNelem];
+  for(Int_t i = 0; i < fNelem; i++) {
+    fPedSum[i]   = 0;
+    fPedSum2[i]  = 0;
     fPedCount[i] = 0;
   }
-
 }
 
 //_____________________________________________________________________________
@@ -481,7 +507,7 @@ void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
     Int_t element = hit->fCounter - 1;
     Int_t nadc = hit->GetRawAdcHitPos().GetPulseIntRaw();
     if(nadc <= fPedLimit[element]) {
-      fPedSum[element] += nadc;
+      fPedSum[element]  += nadc;
       fPedSum2[element] += nadc*nadc;
       fPedCount[element]++;
       if(fPedCount[element] == fMinPeds/5) {
@@ -497,7 +523,7 @@ void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits)
 }
 
 //_____________________________________________________________________________
-void THcCherenkov::CalculatePedestals( )
+void THcCherenkov::CalculatePedestals()
 {
   // Use the accumulated pedestal data to calculate pedestals
   // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
@@ -523,34 +549,30 @@ void THcCherenkov::CalculatePedestals( )
 }
 
 //_____________________________________________________________________________
-Int_t THcCherenkov::GetCerIndex( Int_t nRegion, Int_t nValue ) {
-
-  return fCerNRegions * nValue + nRegion;
+Int_t THcCherenkov::GetIndex(Int_t nRegion, Int_t nValue)
+{
+  return fNRegions * nValue + nRegion;
 }
 
 
 //_____________________________________________________________________________
-void THcCherenkov::Print( const Option_t* opt) const {
+void THcCherenkov::Print(const Option_t* opt) const
+{
   THaNonTrackingDetector::Print(opt);
-
   // Print out the pedestals
-
   cout << endl;
   cout << "Cherenkov Pedestals" << endl;
-
   // Ahmed
   cout << "No.   ADC" << endl;
-  for(Int_t i=0; i<fNelem; i++){
+  for(Int_t i=0; i<fNelem; i++)
     cout << " " << i << "    " << fPed[i] << endl;
-  }
-
   cout << endl;
 }
 
 //_____________________________________________________________________________
-Double_t THcCherenkov::GetCerNPE() {
-
-  return fNPEsum;
+Double_t THcCherenkov::GetCerNPE()
+{
+  return fNpeSum;
 }
 
 ClassImp(THcCherenkov)
diff --git a/src/THcCherenkov.h b/src/THcCherenkov.h
index 2a02df8..80d18e1 100644
--- a/src/THcCherenkov.h
+++ b/src/THcCherenkov.h
@@ -15,85 +15,91 @@
 class THcCherenkov : public THaNonTrackingDetector, public THcHitList {
 
  public:
-  THcCherenkov( const char* name, const char* description = "",
-		THaApparatus* a = NULL );
+  THcCherenkov(const char* name, const char* description = "", THaApparatus* a = NULL);
   virtual ~THcCherenkov();
 
-  virtual void 	     Clear( Option_t* opt="" );
-  virtual Int_t      Decode( const THaEvData& );
-  virtual EStatus    Init( const TDatime& run_time );
-  void               InitArrays();
-  void               DeleteArrays();
-  virtual Int_t      ReadDatabase( const TDatime& date );
-  virtual Int_t      DefineVariables( EMode mode = kDefine );
-  virtual Int_t      CoarseProcess( TClonesArray& tracks );
-  virtual Int_t      FineProcess( TClonesArray& tracks );
-
-  virtual void AccumulatePedestals(TClonesArray* rawhits);
-  virtual void CalculatePedestals();
-
-  virtual Int_t      ApplyCorrections( void );
-
-  virtual void Print(const Option_t* opt) const;
-
-  Int_t GetCerIndex(Int_t nRegion, Int_t nValue);
+  virtual void 	  Clear(Option_t* opt="");
+  virtual void    Print(const Option_t* opt) const;
+  virtual void    AccumulatePedestals(TClonesArray* rawhits);
+  virtual void    CalculatePedestals();
+  virtual Int_t   Decode(const THaEvData&);
+  virtual Int_t   ReadDatabase(const TDatime& date);
+  virtual Int_t   DefineVariables(EMode mode = kDefine);
+  virtual Int_t   CoarseProcess(TClonesArray& tracks);
+  virtual Int_t   FineProcess(TClonesArray& tracks);
+  virtual Int_t   ApplyCorrections( void );
+  virtual EStatus Init(const TDatime& run_time);
+
+  void  InitArrays();
+  void  DeleteArrays();
+  Int_t GetIndex(Int_t nRegion, Int_t nValue);
 
   //  Double_t GetCerNPE() { return fNPEsum;}
   Double_t GetCerNPE();
 
+  // Vector/TClonesArray length parameters
+  static const Int_t MaxNumCerPmt   = 4;
+  static const Int_t MaxNumAdcPulse = 4;
+
   THcCherenkov();  // for ROOT I/O
  protected:
-  Int_t         fAnalyzePedestals;
-
-  // Parameters
-  Double_t*     fGain;
-  Double_t*     fCerWidth;
-
-  // Event information
-  Int_t         fNhits;
-  Int_t*        fADC_hit;         // [fNelem] Array of flag if ADC hit 1 means  
-  Int_t*        fNPMT;            // [fNelem] Array of ADC amplitudes
-  Double_t*     fADC;             // [fNelem] Array of ADC amplitudes
-  Double_t*     fADC_P;           // [fNelem] Array of ADC amplitudes
-  Double_t*     fNPE;             // [fNelem] Array of ADC amplitudes
-  Double_t      fNPEsum;
-  Int_t         fNCherHit;
-
-  Double_t*        fCerRegionValue;
-  Double_t         fCerChi2Max;
-  Double_t         fCerBetaMin;
-  Double_t         fCerBetaMax;
-  Double_t         fCerETMin;
-  Double_t         fCerETMax;
-  Double_t         fCerMirrorZPos;
-  Int_t            fCerNRegions;
-  Int_t            fCerRegionsValueMax;
-  Int_t*           fCerTrackCounter;     // [fCerNRegions] Array of Cher regions
-  Int_t*           fCerFiredCounter;     // [fCerNRegions] Array of Cher regions
-  Double_t         fCerThresh;
-
-  // Hits
-  TClonesArray* fADCHits;
-
-  // Pedestals
-  Int_t         fNPedestalEvents;
-  Int_t         fMinPeds;
-  Int_t*        fPedSum;	  /* Accumulators for pedestals */
-  Int_t*        fPedSum2;
-  Int_t*        fPedLimit;
-  Double_t*     fPedMean; 	  /* Can be supplied in parameters and then */
-  Int_t*        fPedCount;
-  Double_t*     fPed;
-  Double_t*     fThresh;
-
+  Int_t     fAnalyzePedestals;
+  Int_t     fDebugAdc;
+  Double_t* fWidth;
+
+  Int_t     fNhits;
+  Int_t     fTotNumAdcHits;
+  Int_t     fTotNumGoodAdcHits;
+  Int_t     fTotNumTracksMatched;
+  Int_t     fTotNumTracksFired;
+  Double_t  fNpeSum;
+  Double_t* fGain;
+
+  vector<Int_t>    fNumAdcHits;
+  vector<Int_t>    fNumGoodAdcHits;
+  vector<Int_t>    fNumTracksMatched;
+  vector<Int_t>    fNumTracksFired;
+  vector<Double_t> fGoodAdcPed;
+  vector<Double_t> fGoodAdcPulseInt;
+  vector<Double_t> fGoodAdcPulseIntRaw;
+  vector<Double_t> fGoodAdcPulseAmp;
+  vector<Double_t> fGoodAdcPulseTime;
+  vector<Double_t> fNpe;
+
+  Int_t     fNRegions;
+  Int_t     fRegionsValueMax;
+  Double_t  fRedChi2Min;
+  Double_t  fRedChi2Max;
+  Double_t  fBetaMin;
+  Double_t  fBetaMax;
+  Double_t  fENormMin;
+  Double_t  fENormMax;
+  Double_t  fMirrorZPos;
+  Double_t  fNpeThresh;
+  Double_t  fAdcTimeWindowMin;
+  Double_t  fAdcTimeWindowMax;
+  Double_t* fRegionValue;
+
+  // 6 Gev pedestal variables
+  Int_t     fNPedestalEvents;
+  Int_t     fMinPeds;
+  Int_t*    fPedSum;	  /* Accumulators for pedestals */
+  Int_t*    fPedSum2;
+  Int_t*    fPedLimit;
+  Int_t*    fPedCount;
+  Double_t* fPedMean; 	  /* Can be supplied in parameters and then */
+  Double_t* fPed;
+  Double_t* fThresh;
+
+  // 12 Gev FADC variables
   TClonesArray* frAdcPedRaw;
   TClonesArray* frAdcPulseIntRaw;
   TClonesArray* frAdcPulseAmpRaw;
   TClonesArray* frAdcPulseTimeRaw;
-
   TClonesArray* frAdcPed;
   TClonesArray* frAdcPulseInt;
   TClonesArray* frAdcPulseAmp;
+  TClonesArray* fAdcErrorFlag;
 
   void Setup(const char* name, const char* description);
   virtual void  InitializePedestals( );
-- 
GitLab