containers/docker/README.md

-Hall A/C Software
-=================
-
-
-Images for running containers for all aspects of the SANE experimental 
-analysis.
-
-The starting point is a pre-built image for the ROOT libraries. (ubuntu + ROOT)
-
-Main software libraries:
-
- - `evio`: Built from https://github.com/whit2333/hallac_evio
- - `analyzer`: Hall A analyzer (podd)from https://github.com/whit2333/analyzer
- - `hcana`: Hall C analyzer from  https://github.com/whit2333/hcana
-
-These are all built using the super build project `cool_halls` (https://github.com/whit2333/cool_halls)
-
-

containers/docker/deploy.env

-# You have to define the values in {}
-#DOCKER_REPO={account-nr}.dkr.ecr.{region}.amazonaws.com
-## optional aws-cli options
-#AWS_CLI_PROFILE={aws-cli-profile}
-#AWS_CLI_REGION={aws-cli-region}

containers/docker/usage.sh

-# INSTALL
-# - copy the files deploy.env, config.env, version.sh and Makefile to your repo
-# - replace the vars in deploy.env
-# - define the version script
-
-# Build the container
-make build
-
-# Build and publish the container
-make release
-
-# Publish a container to AWS-ECR.
-# This includes the login to the repo
-make publish
-
-# Run the container
-make run
-
-# Build an run the container
-make up
-
-# Stop the running container
-make stop
-
-# Build the container with differnt config and deploy file
-make cnf=another_config.env dpl=another_deploy.env build
\ No newline at end of file

containers/singularity/Singularity

-Bootstrap: docker
-From: eicweb.phy.anl.gov:4567/jlab/hallc/analyzer_software/hcana:latest
-
-%help
-  Hall A/C container.
-  Tools:
-     - evio        : EVIO DAQ data format  
-     - analyzer    : Hall A analyzer (podd) 
-     - hcana       : Hall C analyzer (hcana)
-     - root        : root version used for the analyzer
-     - rootls, rootbrowse, root_config
-
-%labels
-  Maintainer "Whitney Armstrong, Sylvester Joosten"
-  Version v1.0
-
-%setup -c /bin/bash
-  export SINGULARITY_SHELL=/bin/bash
-
-%environment -c /bin/bash
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include:/usr/local/include/podd:/usr/local/include/hcana
-
-%post -c /bin/bash
-  echo "Hello from post"
-  echo "Install additional software here"
-  source /usr/local/bin/thisroot.sh
-  ## libformat and nlohmann json used heavily in new generation replay scripts
-  ## libformat
-  #cd /tmp && git clone https://github.com/fmtlib/fmt.git && cd fmt && \
-  #  git checkout 5.3.0 && mkdir /tmp/build && cd /tmp/build && \
-  #  cmake -DBUILD_SHARED_LIBS=TRUE ../fmt &&
-  #  make -j4 install && cd /tmp && rm -r /tmp/build && rm -r /tmp/fmt
-  ### json
-
-# =======================
-# global
-# =======================
-
-%runscript
-  echo "Launching a shell in the Hall A/C singularity container
-  exec bash
-
-
-# =======================
-# root
-# =======================
-%apprun root
-  root "$@"
-
-%appenv root
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
-
-# =======================
-# analyzer
-# =======================
-%apprun analyzer
-  analyzer "$@"
-
-%appenv analyzer
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
-
-# =======================
-# hcana
-# =======================
-%apphelp hcana
-  Run the Hall-C analyzer with same root-style arguments.
-
-%apprun hcana
-  source /usr/local/bin/thisroot.sh
-  hcana "$@"
-
-%appenv hcana
-  export DB_DIR=DBASE
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOTSYS=/usr/local
-  export ROOT_INCLUDE_PATH=/usr/local/include
-  export ROOT_INCLUDE_PATH=/usr/local/include:/usr/local/include/podd:/usr/local/include/hcana
-
-# =======================
-# root-config
-# =======================
-%apprun root_config
-  root-config "$@"
-
-%appenv root_config
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
-
-
-# =======================
-# rootbrowse
-# =======================
-%apprun rootbrowse
-  rootbrowse "$@"
-
-%appenv rootbrowse
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana
-
-# =======================
-# rootls
-# =======================
-%apprun rootls
-  rootls "$@"
-
-%appenv rootls
-  export PYTHONPATH=/usr/local/lib:${PYTHONPATH}
-  export PATH=/usr/local/bin:${PATH}
-  export LD_LIBRARY_PATH=/usr/local/lib:$LD_LIBRARY_PATH
-  export ROOT_INCLUDE_PATH=/usr/local/include/podd:/usr/local/include/hcana

podd @ 355b3ce6

-Subproject commit ecff1455f2868eee65295d4aa6e772eefb2dd941
+Subproject commit 355b3ce629f8f47f107e2a0fcfa7ee55b3463b6c

src/THcAerogel.cxx

@@ -39,7 +39,8 @@ THcAerogel::THcAerogel( const char* name, const char* description,
   THaNonTrackingDetector(name,description,apparatus),
   fPresentP(0),
   fAdcPosTimeWindowMin(0), fAdcPosTimeWindowMax(0), fAdcNegTimeWindowMin(0),
-  fAdcNegTimeWindowMax(0), fRegionValue(0), fPosGain(0), fNegGain(0),
+  fAdcNegTimeWindowMax(0),fPedNegDefault(0),fPedPosDefault(0),
+  fRegionValue(0), fPosGain(0), fNegGain(0),
   frPosAdcPedRaw(0), frPosAdcPulseIntRaw(0), frPosAdcPulseAmpRaw(0),
   frPosAdcPulseTimeRaw(0), frPosAdcPed(0), frPosAdcPulseInt(0),
   frPosAdcPulseAmp(0), frPosAdcPulseTime(0), frNegAdcPedRaw(0),
@@ -59,7 +60,8 @@ THcAerogel::THcAerogel( const char* name, const char* description,
 THcAerogel::THcAerogel( ) :
   THaNonTrackingDetector(),
   fAdcPosTimeWindowMin(0), fAdcPosTimeWindowMax(0), fAdcNegTimeWindowMin(0),
-  fAdcNegTimeWindowMax(0), fRegionValue(0), fPosGain(0), fNegGain(0),
+  fAdcNegTimeWindowMax(0), 
+  fPedNegDefault(0),fPedPosDefault(0),fRegionValue(0), fPosGain(0), fNegGain(0),
   frPosAdcPedRaw(0), frPosAdcPulseIntRaw(0), frPosAdcPulseAmpRaw(0),
   frPosAdcPulseTimeRaw(0), frPosAdcPed(0), frPosAdcPulseInt(0),
   frPosAdcPulseAmp(0), frPosAdcPulseTime(0), frNegAdcPedRaw(0),
@@ -111,6 +113,8 @@ void THcAerogel::DeleteArrays()
   delete [] fAdcPosTimeWindowMax; fAdcPosTimeWindowMax = 0;
   delete [] fAdcNegTimeWindowMin; fAdcNegTimeWindowMin = 0;
   delete [] fAdcNegTimeWindowMax; fAdcNegTimeWindowMax = 0;
+  delete [] fPedNegDefault; fPedNegDefault = 0;
+  delete [] fPedPosDefault; fPedPosDefault = 0;
 
   // 6 GeV variables
   delete fPosTDCHits; fPosTDCHits = NULL;
@@ -296,6 +300,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
   fAdcPosTimeWindowMax = new Double_t [fNelem];
   fAdcNegTimeWindowMin = new Double_t [fNelem];
   fAdcNegTimeWindowMax = new Double_t [fNelem];
+  fPedNegDefault = new Int_t [fNelem];
+  fPedPosDefault = new Int_t [fNelem];
 
   DBRequest list[]={
     {"aero_num_regions",      &fNRegions,         kInt},
@@ -315,6 +321,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
     {"aero_adcPosTimeWindowMax", fAdcPosTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem), 1},
     {"aero_adcNegTimeWindowMin", fAdcNegTimeWindowMin, kDouble, static_cast<UInt_t>(fNelem), 1},
     {"aero_adcNegTimeWindowMax", fAdcNegTimeWindowMax, kDouble, static_cast<UInt_t>(fNelem), 1},
+    {"aero_PedNegDefault", fPedNegDefault, kInt, static_cast<UInt_t>(fNelem), 1},
+    {"aero_PedPosDefault", fPedPosDefault, kInt, static_cast<UInt_t>(fNelem), 1},
     {"aero_adc_tdc_offset",   &fAdcTdcOffset,     kDouble, 0, 1},
     {"aero_debug_adc",        &fDebugAdc,         kInt,    0, 1},
     {"aero_six_gev_data",     &fSixGevData,       kInt,    0, 1},
@@ -331,6 +339,8 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
    fAdcNegTimeWindowMin[ip] = -1000.;
    fAdcPosTimeWindowMax[ip] = 1000.;
    fAdcNegTimeWindowMax[ip] = 1000.;
+   fPedNegDefault[ip] = 0.;
+   fPedPosDefault[ip] = 0.;
   }
 
   fSixGevData = 0; // Set 6 GeV data parameter to false unless set in parameter file
@@ -646,6 +656,21 @@ Int_t THcAerogel::Decode( const THaEvData& evdata )
 
       if (rawPosAdcHit.GetPulseAmpRaw(thit) > 0)  ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 0);
       if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(nrPosAdcHits))->Set(npmt, 1);
+ 
+     if (rawPosAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawPosAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawPosAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawPosAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawPosAdcHit.GetAdcTomV();
+	if (fPedPosDefault[npmt-1] !=0) {
+	  Double_t tPulseInt = AdcToC*(rawPosAdcHit.GetPulseIntRaw(thit) - fPedPosDefault[npmt-1]*PeakPedRatio);
+	  ((THcSignalHit*) frPosAdcPulseInt->ConstructedAt(nrPosAdcHits))->Set(npmt, tPulseInt);
+          ((THcSignalHit*) frPosAdcPedRaw->ConstructedAt(nrPosAdcHits))->Set(npmt, fPedPosDefault[npmt-1]);
+          ((THcSignalHit*) frPosAdcPed->ConstructedAt(nrPosAdcHits))->Set(npmt, float(fPedPosDefault[npmt-1])/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(nrPosAdcHits))->Set(npmt, 0.);	
+      }
 
       ++nrPosAdcHits;
       fTotNumAdcHits++;
@@ -669,6 +694,21 @@ Int_t THcAerogel::Decode( const THaEvData& evdata )
       if (rawNegAdcHit.GetPulseAmpRaw(thit) > 0)  ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 0);
       if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(nrNegAdcHits))->Set(npmt, 1);
 
+     if (rawNegAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawNegAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawNegAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawNegAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawNegAdcHit.GetAdcTomV();
+	if (fPedNegDefault[npmt-1] !=0) {
+	  Double_t tPulseInt = AdcToC*(rawNegAdcHit.GetPulseIntRaw(thit) - fPedNegDefault[npmt-1]*PeakPedRatio);
+	  ((THcSignalHit*) frNegAdcPulseInt->ConstructedAt(nrNegAdcHits))->Set(npmt, tPulseInt);
+          ((THcSignalHit*) frNegAdcPedRaw->ConstructedAt(nrNegAdcHits))->Set(npmt, fPedNegDefault[npmt-1]);
+          ((THcSignalHit*) frNegAdcPed->ConstructedAt(nrNegAdcHits))->Set(npmt, float(fPedNegDefault[npmt-1])/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(nrNegAdcHits))->Set(npmt, 0.);	
+      }
+
       ++nrNegAdcHits;
       fTotNumAdcHits++;
       fTotNumNegAdcHits++;
@@ -690,6 +730,8 @@ Int_t THcAerogel::CoarseProcess( TClonesArray&  ) //tracks
 {
   Double_t StartTime = 0.0;
   if( fglHod ) StartTime = fglHod->GetStartTime();
+   Double_t OffsetTime = 0.0;
+   if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
   //cout << " starttime = " << StartTime << endl;
     // Loop over the elements in the TClonesArray
     for(Int_t ielem = 0; ielem < frPosAdcPulseInt->GetEntries(); ielem++) {
@@ -700,18 +742,13 @@ Int_t THcAerogel::CoarseProcess( TClonesArray&  ) //tracks
       Double_t pulseIntRaw  = ((THcSignalHit*) frPosAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
       Double_t pulseAmp     = ((THcSignalHit*) frPosAdcPulseAmp->ConstructedAt(ielem))->GetData();
       Double_t pulseTime    = ((THcSignalHit*) frPosAdcPulseTime->ConstructedAt(ielem))->GetData();
-      Double_t adctdcdiffTime = StartTime-pulseTime;
-      Bool_t   errorFlag    = ((THcSignalHit*) fPosAdcErrorFlag->ConstructedAt(ielem))->GetData();
+      Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
       ////      Bool_t   pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
       Bool_t   pulseTimeCut = adctdcdiffTime > fAdcPosTimeWindowMin[npmt] && adctdcdiffTime < fAdcPosTimeWindowMax[npmt];
 
       // By default, the last hit within the timing cut will be considered "good"
-     if (!errorFlag)
-      {
 	fGoodPosAdcMult.at(npmt) += 1;
-      }
-
-     if (!errorFlag && pulseTimeCut) {
+     if (pulseTimeCut) {
     	fGoodPosAdcPed.at(npmt)         = pulsePed;
  	//	cout << " out = " << npmt << " " <<   frPosAdcPulseInt->GetEntries() << " " <<fGoodPosAdcMult.at(npmt); 
     	fGoodPosAdcPulseInt.at(npmt)    = pulseInt;
@@ -738,17 +775,13 @@ Int_t THcAerogel::CoarseProcess( TClonesArray&  ) //tracks
       Double_t pulseIntRaw  = ((THcSignalHit*) frNegAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
       Double_t pulseAmp     = ((THcSignalHit*) frNegAdcPulseAmp->ConstructedAt(ielem))->GetData();
       Double_t pulseTime    = ((THcSignalHit*) frNegAdcPulseTime->ConstructedAt(ielem))->GetData();
-      Double_t adctdcdiffTime = StartTime-pulseTime;
-      Bool_t   errorFlag    = ((THcSignalHit*) fNegAdcErrorFlag->ConstructedAt(ielem))->GetData();
+      Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
       ////      Bool_t   pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin && adctdcdiffTime < fAdcTimeWindowMax;
       Bool_t   pulseTimeCut = adctdcdiffTime > fAdcNegTimeWindowMin[npmt] && adctdcdiffTime < fAdcNegTimeWindowMax[npmt];
-      if (!errorFlag)
-      {
 	fGoodNegAdcMult.at(npmt) += 1;
-      }
  
       // By default, the last hit within the timing cut will be considered "good"
-      if (!errorFlag && pulseTimeCut) {
+      if (pulseTimeCut) {
     	fGoodNegAdcPed.at(npmt)         = pulsePed;
     	fGoodNegAdcPulseIntRaw.at(npmt) = pulseIntRaw;
     	fGoodNegAdcPulseAmp.at(npmt)    = pulseAmp;

src/THcAerogel.h

@@ -73,6 +73,8 @@ class THcAerogel : public THaNonTrackingDetector, public THcHitList {
   Double_t  *fAdcPosTimeWindowMax;
   Double_t  *fAdcNegTimeWindowMin;
   Double_t  *fAdcNegTimeWindowMax;
+  Int_t*  fPedNegDefault;
+  Int_t*  fPedPosDefault;
 
   Double_t  fAdcTdcOffset;
   Double_t  *fRegionValue;

src/THcCherenkov.cxx

@@ -129,27 +129,20 @@ void THcCherenkov::DeleteArrays() {
   fGain = NULL;
 
   // 6 Gev variables
-  delete[] fPedSum;
-  fPedSum = NULL;
-  delete[] fPedSum2;
-  fPedSum2 = NULL;
-  delete[] fPedLimit;
-  fPedLimit = NULL;
-  delete[] fPedMean;
-  fPedMean = NULL;
-  delete[] fPedCount;
-  fPedCount = NULL;
-  delete[] fPed;
-  fPed = NULL;
-  delete[] fThresh;
-  fThresh = NULL;
-
-  delete[] fAdcTimeWindowMin;
-  fAdcTimeWindowMin = 0;
-  delete[] fAdcTimeWindowMax;
-  fAdcTimeWindowMax = 0;
-  delete[] fRegionValue;
-  fRegionValue = 0;
+
+  delete [] fPedSum;   fPedSum   = NULL;
+  delete [] fPedSum2;  fPedSum2  = NULL;
+  delete [] fPedLimit; fPedLimit = NULL;
+  delete [] fPedMean;  fPedMean  = NULL;
+  delete [] fPedCount; fPedCount = NULL;
+  delete [] fPed;      fPed      = NULL;
+  delete [] fThresh;   fThresh   = NULL;
+
+  delete [] fPedDefault; fPedDefault = 0;
+  delete [] fAdcTimeWindowMin; fAdcTimeWindowMin = 0;
+  delete [] fAdcTimeWindowMax; fAdcTimeWindowMax = 0;
+  delete [] fRegionValue; fRegionValue = 0;
+
 }
 
 //_____________________________________________________________________________
@@ -212,39 +205,48 @@ Int_t THcCherenkov::ReadDatabase(const TDatime& date) {
   //     << GetName() << " with " << fNelem << " PMTs" << endl;
 
   // 6 GeV pedestal paramters
-  fPedLimit         = new Int_t[fNelem];
-  fGain             = new Double_t[fNelem];
-  fPedMean          = new Double_t[fNelem];
-  fAdcTimeWindowMin = new Double_t[fNelem];
-  fAdcTimeWindowMax = new Double_t[fNelem];
+
+  fPedLimit = new Int_t[fNelem];
+  fGain     = new Double_t[fNelem];
+  fPedMean  = new Double_t[fNelem];
+  fAdcTimeWindowMin =  new Double_t[fNelem];
+  fAdcTimeWindowMax=  new Double_t[fNelem];
+  fPedDefault=  new Int_t[fNelem];
+
   // Region parameters
   fRegionsValueMax = fNRegions * 8;
   fRegionValue     = new Double_t[fRegionsValueMax];
   fAdcGoodElem = new Int_t[fNelem];
   fAdcPulseAmpTest = new Double_t[fNelem];
 
-  DBRequest list[] = {{"_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},
-                      {"_dp_min", &fDpMin, kDouble},
-                      {"_dp_max", &fDpMax, kDouble},
-                      {"_mirror_zpos", &fMirrorZPos, kDouble},
-                      {"_npe_thresh", &fNpeThresh, kDouble},
-                      {"_debug_adc", &fDebugAdc, kInt, 0, 1},
-                      {"_adcTimeWindowMin", fAdcTimeWindowMin, kDouble, (UInt_t)fNelem, 1},
-                      {"_adcTimeWindowMax", fAdcTimeWindowMax, kDouble, (UInt_t)fNelem, 1},
-                      {"_adc_tdc_offset", &fAdcTdcOffset, kDouble, 0, 1},
-                      {"_region", &fRegionValue[0], kDouble, (UInt_t)fRegionsValueMax},
-                      {"_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
-                      {0}};
-  for (Int_t i = 0; i < fNelem; i++) {
-    fAdcTimeWindowMin[i] = -1000.;
-    fAdcTimeWindowMax[i] = 1000.;
+
+  DBRequest list[]={
+    {"_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},
+    {"_dp_min",           &fDpMin,            kDouble},
+    {"_dp_max",           &fDpMax,            kDouble},
+    {"_mirror_zpos",      &fMirrorZPos,       kDouble},
+    {"_npe_thresh",       &fNpeThresh,        kDouble},
+    {"_debug_adc",        &fDebugAdc,         kInt, 0, 1},
+    {"_adcTimeWindowMin", fAdcTimeWindowMin, kDouble,(UInt_t) fNelem,1},
+    {"_adcTimeWindowMax", fAdcTimeWindowMax, kDouble, (UInt_t) fNelem,1},
+    {"_PedDefault", fPedDefault, kInt, (UInt_t) fNelem,1},
+    {"_adc_tdc_offset",   &fAdcTdcOffset,     kDouble, 0, 1},
+    {"_region",           &fRegionValue[0],   kDouble,  (UInt_t) fRegionsValueMax},
+    {"_adcrefcut",        &fADC_RefTimeCut,   kInt,    0, 1},
+    {0}
+  };
+  for (Int_t i=0;i<fNelem;i++) {
+    fAdcTimeWindowMin[i]=-1000.;
+    fAdcTimeWindowMax[i]=1000.;
+    fPedDefault[i]=0;
+
   }
   fDebugAdc       = 0; // Set ADC debug parameter to false unless set in parameter file
   fAdcTdcOffset   = 0.0;
@@ -298,32 +300,36 @@ Int_t THcCherenkov::DefineVariables(EMode mode) {
   } // end debug statement
 
   RVarDef vars[] = {
-      {"adcCounter", "ADC counter numbers", "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
-      {"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"},
-      {"numGoodAdcHits", "Number of Good ADC Hits Per PMT",
-       "fNumGoodAdcHits"}, // Cherenkov occupancy
-      {"totNumGoodAdcHits", "Total Number of Good ADC Hits",
-       "fTotNumGoodAdcHits"}, // Cherenkov multiplicity
-      {"numTracksMatched", "Number of Tracks Matched Per Region", "fNumTracksMatched"},
-      {"numTracksFired", "Number of Tracks that Fired Per Region", "fNumTracksFired"},
-      {"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"},
-      {"totNumTracksFired", "Total Number of Tracks that Fired", "fTotNumTracksFired"},
-
-      {"xAtCer", "Track X at Cherenkov mirror", "fXAtCer"},
-      {"yAtCer", "Track Y at Cherenkov mirror", "fYAtCer"},
-
-      {"npe", "Number of PEs", "fNpe"},
-      {"npeSum", "Total Number of PEs", "fNpeSum"},
-
-      {"goodAdcPed", "Good ADC pedestals", "fGoodAdcPed"},
-      {"goodAdcMult", "Good ADC Multiplicity", "fGoodAdcMult"},
-      {"goodAdcHitUsed", "Good ADC Hit Used", "fGoodAdcHitUsed"},
-      {"goodAdcPulseInt", "Good ADC pulse integrals", "fGoodAdcPulseInt"},
-      {"goodAdcPulseIntRaw", "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"},
-      {"goodAdcPulseAmp", "Good ADC pulse amplitudes", "fGoodAdcPulseAmp"},
-      {"goodAdcPulseTime", "Good ADC pulse times", "fGoodAdcPulseTime"},
-      {"goodAdcTdcDiffTime", "Good Hodo Start - ADC pulse times", "fGoodAdcTdcDiffTime"},
-      {0}};
+
+    {"adcCounter",   "ADC counter numbers",            "frAdcPulseIntRaw.THcSignalHit.GetPaddleNumber()"},
+    {"adcErrorFlag", "Error Flag for When FPGA Fails", "fAdcErrorFlag.THcSignalHit.GetData()"},
+
+    {"numGoodAdcHits",    "Number of Good ADC Hits Per PMT", "fNumGoodAdcHits"},    // Cherenkov occupancy
+    {"totNumGoodAdcHits", "Total Number of Good ADC Hits",   "fTotNumGoodAdcHits"}, // Cherenkov multiplicity
+
+    {"numTracksMatched",    "Number of Tracks Matched Per Region",       "fNumTracksMatched"},
+    {"numTracksFired",      "Number of Tracks that Fired Per Region",    "fNumTracksFired"},
+    {"totNumTracksMatched", "Total Number of Tracks Matched Per Region", "fTotNumTracksMatched"},
+    {"totNumTracksFired",   "Total Number of Tracks that Fired",         "fTotNumTracksFired"},
+
+    {"xAtCer",       "Track X at Cherenkov mirror",    "fXAtCer"},
+    {"yAtCer",       "Track Y at Cherenkov mirror",    "fYAtCer"},
+
+    {"npe",          "Number of PEs",                  "fNpe"},
+    {"npeSum",       "Total Number of PEs",            "fNpeSum"},
+
+    {"goodAdcPed",          "Good ADC pedestals",           "fGoodAdcPed"},
+    {"goodAdcMult",          "Good ADC Multiplicity",           "fGoodAdcMult"},
+    {"goodAdcHitUsed",          "Good ADC Hit Used",           "fGoodAdcHitUsed"},
+    {"goodAdcPulseInt",     "Good ADC pulse integrals",     "fGoodAdcPulseInt"},
+    {"goodAdcPulseIntRaw",  "Good ADC raw pulse integrals", "fGoodAdcPulseIntRaw"},
+    {"goodAdcPulseAmp",     "Good ADC pulse amplitudes",    "fGoodAdcPulseAmp"},
+    {"goodAdcPulseTime",    "Good ADC pulse times",         "fGoodAdcPulseTime"},
+     {"goodAdcTdcDiffTime",    "Good Hodo Start - ADC pulse times",         "fGoodAdcTdcDiffTime"},
+       {"RefTime",      "Raw ADC RefTime (chan) ", "fRefTime"},        // Raw reference time
+  { 0 }
+  };
+
 
   return DefineVarsFromList(vars, mode);
 }
@@ -360,6 +366,7 @@ void THcCherenkov::Clear(Option_t* opt) {
   fYAtCer = 0.0;
 
   fNpeSum = 0.0;
+  fRefTime=kBig;
 
   frAdcPedRaw->Clear();
   frAdcPulseIntRaw->Clear();
@@ -406,7 +413,12 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
   }
   fNhits = DecodeToHitList(evdata, !present);
 
-  if (gHaCuts->Result("Pedestal_event")) {
+
+  //THcHallCSpectrometer *app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
+  // cout << "Cerenkov  Event num = " << evdata.GetEvNum() << " spec = " << app->GetName() << endl;
+
+  if(gHaCuts->Result("Pedestal_event")) {
+
     AccumulatePedestals(fRawHitList);
     fAnalyzePedestals = 1; // Analyze pedestals first normal events
     return (0);
@@ -421,14 +433,16 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
   UInt_t nrAdcHits = 0;
   _waveforms.clear();
 
-  while (ihit < fNhits) {
-
-    THcCherenkovHit* hit       = (THcCherenkovHit*)fRawHitList->At(ihit);
-    Int_t            npmt      = hit->fCounter;
-    THcRawAdcHit&    rawAdcHit = hit->GetRawAdcHitPos();
-    
-    _waveforms.push_back({rawAdcHit.GetSampleBuffer()});
 
+  while(ihit < fNhits) {
+ 
+    THcCherenkovHit* hit         = (THcCherenkovHit*) fRawHitList->At(ihit);
+    Int_t            npmt        = hit->fCounter;
+    THcRawAdcHit&    rawAdcHit   = hit->GetRawAdcHitPos();
+    if (rawAdcHit.GetNPulses() >0 && rawAdcHit.HasRefTime()) {
+      fRefTime=rawAdcHit.GetRefTime() ;
+    }
+    //if (rawAdcHit.GetNPulses()>0) cout << "Cer npmt = " << " ped = " << rawAdcHit.GetPed() << endl;
 
     for (UInt_t thit = 0; thit < rawAdcHit.GetNPulses(); thit++) {
 
@@ -455,6 +469,23 @@ Int_t THcCherenkov::Decode(const THaEvData& evdata) {
       if (rawAdcHit.GetPulseAmpRaw(thit) <= 0)
         ((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 1);
 
+      if (rawAdcHit.GetPulseAmpRaw(thit) <= 0) {
+	Double_t PeakPedRatio= rawAdcHit.GetF250_PeakPedestalRatio();
+	Int_t NPedSamples= rawAdcHit.GetF250_NPedestalSamples();
+	Double_t AdcToC =  rawAdcHit.GetAdcTopC();
+	Double_t AdcToV =  rawAdcHit.GetAdcTomV();
+	if (fPedDefault[npmt-1] !=0) {
+	  Double_t tPulseInt = AdcToC*(rawAdcHit.GetPulseIntRaw(thit) - fPedDefault[npmt-1]*PeakPedRatio);
+	  ((THcSignalHit*) frAdcPulseInt->ConstructedAt(nrAdcHits))->Set(npmt, tPulseInt);
+          ((THcSignalHit*) frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(npmt, fPedDefault[npmt-1]);
+          ((THcSignalHit*) frAdcPed->ConstructedAt(nrAdcHits))->Set(npmt, float(fPedDefault[npmt-1])/float(NPedSamples)*AdcToV);
+	  
+	}
+	((THcSignalHit*) frAdcPulseAmp->ConstructedAt(nrAdcHits))->Set(npmt, 0.);
+	
+      }
+
+	
       ++nrAdcHits;
       fTotNumAdcHits++;
       fNumAdcHits.at(npmt - 1) = npmt;
@@ -471,7 +502,9 @@ Int_t THcCherenkov::ApplyCorrections(void) { return (0); }
 Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
   Double_t StartTime = 0.0;
   if( fglHod ) StartTime = fglHod->GetStartTime();
-  for(Int_t ipmt = 0; ipmt < fNelem; ipmt++) {
+   Double_t OffsetTime = 0.0;
+   if( fglHod ) OffsetTime = fglHod->GetOffsetTime();
+ for(Int_t ipmt = 0; ipmt < fNelem; ipmt++) {
     fAdcPulseAmpTest[ipmt] = -1000.;
     fAdcGoodElem[ipmt]=-1;
    }
@@ -480,17 +513,18 @@ Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
     Int_t    npmt         = ((THcSignalHit*) frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
     Double_t pulseTime    = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
     Double_t pulseAmp     = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
-   Double_t adctdcdiffTime = StartTime-pulseTime;
-     Bool_t   errorFlag    = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
+    Bool_t   errorFlag    = ((THcSignalHit*) fAdcErrorFlag->ConstructedAt(ielem))->GetData();
+    Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
     Bool_t   pulseTimeCut = adctdcdiffTime > fAdcTimeWindowMin[npmt] && adctdcdiffTime < fAdcTimeWindowMax[npmt];
-    if (!errorFlag)
-      {
-	fGoodAdcMult.at(npmt) += 1;
-      }
-    if (!errorFlag && pulseTimeCut && pulseAmp > fAdcPulseAmpTest[npmt]) {
-       fAdcGoodElem[npmt]=ielem;
-       fAdcPulseAmpTest[npmt] = pulseAmp;
-    }
+ 	fGoodAdcMult.at(npmt) += 1;
+	if (!errorFlag) {
+	  if (pulseTimeCut && pulseAmp > fAdcPulseAmpTest[npmt]) {
+             fAdcGoodElem[npmt]=ielem;
+              fAdcPulseAmpTest[npmt] = pulseAmp;
+	  }
+        } else {
+	  if (pulseTimeCut) fAdcGoodElem[npmt]=ielem;
+        }
   }
   // Loop over the npmt
   for(Int_t npmt = 0; npmt < fNelem; npmt++) {
@@ -501,17 +535,18 @@ Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
     Double_t pulseIntRaw  = ((THcSignalHit*) frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
     Double_t pulseAmp     = ((THcSignalHit*) frAdcPulseAmp->ConstructedAt(ielem))->GetData();
     Double_t pulseTime    = ((THcSignalHit*) frAdcPulseTime->ConstructedAt(ielem))->GetData();
-   Double_t adctdcdiffTime = StartTime-pulseTime;
-    // By default, the last hit within the timing cut will be considered "good"
+    Double_t adctdcdiffTime = StartTime-pulseTime+OffsetTime;
       fGoodAdcPed.at(npmt)         = pulsePed;
       fGoodAdcHitUsed.at(npmt)     = ielem + 1;
       fGoodAdcPulseInt.at(npmt)    = pulseInt;
       fGoodAdcPulseIntRaw.at(npmt) = pulseIntRaw;
       fGoodAdcPulseAmp.at(npmt)    = pulseAmp;
       fGoodAdcPulseTime.at(npmt)   = pulseTime;
-      fGoodAdcTdcDiffTime.at(npmt) = adctdcdiffTime;
 
-      fNpe.at(npmt) = fGain[npmt] * fGoodAdcPulseInt.at(npmt);
+      fGoodAdcTdcDiffTime.at(npmt)   = adctdcdiffTime;
+      
+      fNpe.at(npmt) = fGain[npmt]*pulseInt;
+
       fNpeSum += fNpe.at(npmt);
 
       fTotNumGoodAdcHits++;

src/THcCherenkov.h

@@ -69,6 +69,7 @@ public:
   Int_t     fTotNumGoodAdcHits;
   Int_t     fTotNumTracksMatched;
   Int_t     fTotNumTracksFired;
+  Double_t  fRefTime;
   Double_t  fNpeSum;
   Double_t* fGain;
 
@@ -102,6 +103,7 @@ public:
   Double_t  fNpeThresh;
   Double_t*  fAdcTimeWindowMin;
   Double_t*  fAdcTimeWindowMax;
+  Int_t*  fPedDefault;
   Double_t  fAdcTdcOffset;
   Double_t* fRegionValue;
 

src/THcCoinTime.cxx

@@ -65,6 +65,14 @@ void THcCoinTime::Clear( Option_t* opt )
   fROC2_ePosCoinTime=kBig;
   fROC1_RAW_CoinTime=kBig;
   fROC2_RAW_CoinTime=kBig;
+  fTRIG1_ePosCoinTime=kBig;
+  fTRIG4_ePosCoinTime=kBig;
+  fTRIG1_ePiCoinTime=kBig;
+  fTRIG4_ePiCoinTime=kBig;
+  fTRIG1_eKCoinTime=kBig;
+  fTRIG4_eKCoinTime=kBig;
+  fTRIG1_epCoinTime=kBig;
+  fTRIG4_epCoinTime=kBig;
 }
 
 //_____________________________________________________________________________
@@ -263,6 +271,9 @@ Int_t THcCoinTime::Process( const THaEvData& evdata )
       Double_t hms_ypfp = theHMSTrack->GetPhi();        
       Double_t HMS_FPtime = theHMSTrack->GetFPTime();    
       
+      if (SHMS_FPtime==-2000 || HMS_FPtime==-2000)  return 1;
+      if (SHMS_FPtime==-1000 || HMS_FPtime==-1000)  return 1;
+      
       //Get raw TDC Times for HMS/SHMS (3/4 trigger)
       pTRIG1_TdcTime_ROC1 = fCoinDet->Get_CT_Trigtime(0);  //SHMS
       pTRIG4_TdcTime_ROC1 = fCoinDet->Get_CT_Trigtime(1);  //HMS

src/THcDriftChamberPlane.cxx

@@ -298,6 +298,7 @@ Int_t THcDriftChamberPlane::DefineVariables( EMode mode )
     {"dist","Drift distancess",
      "fHits.THcDCHit.GetDist()"},
     {"nhit", "Number of hits", "GetNHits()"},
+    {"RefTime", "TDC reference time", "fTdcRefTime"},
     { 0 }
   };
 
@@ -351,7 +352,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
 
   fHits->Clear();
   fRawHits->Clear();
-
+  fTdcRefTime = kBig;
   Int_t nrawhits = rawhits->GetLast()+1;
   fNRawhits=0;
   Int_t ihit = nexthit;
@@ -365,6 +366,7 @@ Int_t THcDriftChamberPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
     Int_t wireNum = hit->fCounter;
     THcDCWire* wire = GetWire(wireNum);
     Bool_t First_Hit_In_Window = kTRUE;
+    if (hit->GetRawTdcHit().HasRefTime()) fTdcRefTime = hit->GetRawTdcHit().GetRefTime();
     for(UInt_t mhit=0; mhit<hit->GetRawTdcHit().GetNHits(); mhit++) {
       fNRawhits++;
       /* Sort into early, late and ontime */
@@ -391,6 +393,7 @@ Int_t THcDriftChamberPlane::SubtractStartTime()
 {
   Double_t StartTime = 0.0;
   if( fglHod ) StartTime = fglHod->GetStartTime();
+  if (StartTime == -1000) StartTime = 0.0;
   for(Int_t ihit=0;ihit<GetNHits();ihit++) { 
     THcDCHit *thishit = (THcDCHit*) fHits->At(ihit);
     Double_t temptime= thishit->GetTime()-StartTime;

src/THcDriftChamberPlane.h

@@ -86,6 +86,8 @@ protected:
   TClonesArray* fHits;
   TClonesArray* fRawHits;
   TClonesArray* fWires;
+ 
+  Double_t fTdcRefTime;
 
   Int_t fVersion;
   Int_t fWireOrder;

src/THcHallCSpectrometer.cxx

@@ -142,6 +142,10 @@ Int_t THcHallCSpectrometer::DefineVariables( EMode mode )
   fIsSetup = ( mode == kDefine );
   RVarDef vars[] = {
     { "tr.betachisq", "Chi2 of beta", "fTracks.THaTrack.GetBetaChi2()"},
+    { "tr.GoodPlane4", "Flag for track hitting hodo plane 4", "fTracks.THaTrack.GetGoodPlane4()"},
+    { "tr.GoodPlane3", "Flag for track hitting hodo plane 3", "fTracks.THaTrack.GetGoodPlane3()"},
+    { "tr.fptime", "Track hodo focal plane time", "fTracks.THaTrack.GetFPTime()"},
+    { "tr.npmt", "Track number of hodo PMTs hit", "fTracks.THaTrack.GetNPMT()"},
     { "tr.PruneSelect", "Prune Select ID", "fPruneSelect"},
     { "present", "Trigger Type includes this spectrometer", "fPresent"},
     { 0 }

src/THcHelicity.cxx

@@ -32,6 +32,8 @@ THcHelicity::THcHelicity(const char* name, const char* description, THaApparatus
   //  for( Int_t i = 0; i < NHIST; ++i )
   //    fHisto[i] = 0;
   //  memset(fHbits, 0, sizeof(fHbits));
+  fglHelicityScaler = 0;
+  fHelicityHistory  = 0;
 }
 
 //_____________________________________________________________________________
@@ -91,6 +93,8 @@ THaAnalysisObject::EStatus THcHelicity::Init(const TDatime& date) {
   fglHelicityScaler = 0;
   fHelicityHistory  = 0;
 
+  fRecommendedFreq = -1.0;
+
   fStatus = kOK;
   return fStatus;
 }
@@ -120,9 +124,12 @@ Int_t THcHelicity::ReadDatabase(const TDatime& date) {
   fSign                      = 1; // Default helicity sign
   fRingSeed_reported_initial = 0; // Initial see that should predict reported
                                   // helicity of first quartet.
-  fFirstCycle = -1;               // First Cycle that starts a quad (0 to 3)
-  fFreq       = 29.5596;
-  fHelDelay   = 8;
+  fFirstCycle   = -1;             // First Cycle that starts a quad (0 to 3)
+  fNLastQuartet = -1;
+  fNQuartet     = 0;
+  //  fFreq = 29.5596;
+  fFreq     = 120.0007547169;
+  fHelDelay = 8;
 
   DBRequest list[] = {//	  {"_hsign", &fSign, kInt, 0, 1},
                       {"helicity_delay", &fHelDelay, kInt, 0, 1},
@@ -270,42 +277,6 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
   fMPS              = fIsMPS ? 1 : 0;
   fQrt              = fIsQrt ? 1 : 0; // Last of quartet
 
-#if 0
-  if (fglHelicityScaler) {
-    Int_t nhelev  = fglHelicityScaler->GetNevents();
-    Int_t ncycles = fglHelicityScaler->GetNcycles();
-    if (nhelev > 0) {
-      for (Int_t i = 0; i < nhelev; i++) {
-        fScaleQuartet = (fHelicityHistory[i] & 2) != 0;
-        Int_t ispos   = fHelicityHistory[i] & 1;
-        if (fScaleQuartet) {
-          fScalerSeed = ((fScalerSeed << 1) | ispos) & 0x3FFFFFFF;
-          if (fNBits >= fMAXBIT) {
-            Int_t seedscan = fScalerSeed;
-            Int_t nbehind;
-            for (nbehind = 0; nbehind < 4; nbehind++) {
-              if (seedscan == fRingSeed_reported) {
-                if (nbehind > 1) {
-                  cout << "Scaler seed behind " << nbehind << " quartets" << endl;
-                  cout << "Ev seed     " << bitset<32>(fRingSeed_reported) << endl;
-                  cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
-                }
-                break;
-              }
-              seedscan = RanBit30(seedscan);
-            }
-            if (nbehind > 4) {
-              cout << "Scaler seed does not match" << endl;
-              cout << "Ev seed     " << bitset<32>(fRingSeed_reported) << endl;
-              cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
-            }
-          }
-        }
-      }
-    }
-  }
-#endif
-
   if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
     fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
     return 0;
@@ -336,11 +307,17 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
     Int_t missed = 0;
     //    Double_t elapsed_time = (fTITime - fFirstEvTime)/250000000.0;
     if (fIsMPS) {
-      fPeriodCheck       = fmod(fTITime / fTIPeriod, 1.0);
+      fPeriodCheck = fmod(fTITime / (250000000.0 / fFreq) - fPeriodCheckOffset, 1.0);
+
       fCycle             = (fTITime / fTIPeriod);
       fActualHelicity    = kUnknown;
       fPredictedHelicity = kUnknown;
       if (fFoundMPS) {
+        if (fRecommendedFreq < 0.0) {
+          if (TMath::Abs(fPeriodCheck - 0.5) > 0.25) {
+            fRecommendedFreq = fFreq * (1 - (fPeriodCheck - 0.5) / fCycle);
+          }
+        }
         missed = TMath::Nint(fTITime / fTIPeriod - fLastMPSTime / fTIPeriod);
         if (missed < 1) { // was <=1
           fLastMPSTime       = (fTITime + fLastMPSTime + missed * fTIPeriod) / 2;
@@ -355,8 +332,9 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
         //      cout << "Found MPS" << endl;
         // check for Nint((time-last)/period) > 1
       } else {
-        fFoundMPS    = kTRUE;
-        fLastMPSTime = fTITime;
+        fFoundMPS          = kTRUE;
+        fLastMPSTime       = fTITime;
+        fPeriodCheckOffset = (fPeriodCheck - .5);
       }
     } else if (fFoundMPS) {                     //
       if (fTITime - fLastMPSTime > fTIPeriod) { // We missed MPS periods
@@ -367,14 +345,14 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
           Int_t quartets_missed = (newNCycle - fFirstCycle) / 4 - (fNCycle - fFirstCycle) / 4;
           int   quartetphase    = (newNCycle - fFirstCycle) % 4;
           //	  cout << "  " << fNCycle << " " << newNCycle << " " << fFirstCycle << " " <<
-          // quartets_missed << " " << quartetphase << endl; 	  cout << "Cycles " << fNCycle << "
-          // " << newNCycle << " " << fFirstCycle
+          // quartets_missed << " " << quartetphase << endl;
+          //	  cout << "Cycles " << fNCycle << " " << newNCycle << " " << fFirstCycle
           //	       << " skipped " << quartets_missed << " quartets" << endl;
           fNCycle = newNCycle;
           // Need to reset fQuartet to reflect where we are based on the current
           // reported helicity.  So we don't fail quartet testing.
           // But only do this if we are calibrated.
-          if (fNBits >= fMAXBIT) {
+          if (fNBits >= fMAXBIT + 0) {
 
             for (Int_t i = 0; i < quartets_missed; i++) { // Advance the seeds.
               //	      cout << "Advancing seed A " << fNBits << endl;
@@ -407,7 +385,7 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
         fIsNewCycle           = kTRUE;
         fLastReportedHelicity = fReportedHelicity;
       } else { // No missed periods.  Get helicities from rings
-        if (fNBits >= fMAXBIT) {
+        if (fNBits >= fMAXBIT + 0) {
           int quartetphase               = (fNCycle - fFirstCycle) % 4;
           fQuartetStartHelicity          = (fRingSeed_actual & 1) ? kPlus : kMinus;
           fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
@@ -490,17 +468,15 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
                 }
               } else {
                 if (fNCycle - fFirstCycle > 4) { // Not at start of run.  Reset
-                  _logger->warn("Lost quarted sync at cycle {}", fNCycle);
-                  cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " "
-                       << fQuartet[3] << endl;
+                  _logger->warn("Lost quartet sync at cycle {}: {} {} {} {}", fNCycle, fQuartet[0],
+                                fQuartet[1], fQuartet[2], fQuartet[3]);
                   fFirstCycle +=
                       4 * ((fNCycle - fFirstCycle) / 4); // Update, but don't change phase
                 }
                 fFoundQuartet = kFALSE;
                 fNBits        = 0;
-                _logger->info("Searching for first of a quaret at cycle {}", fFirstCycle);
-                cout << fQuartet[0] << " " << fQuartet[1] << " " << fQuartet[2] << " "
-                     << fQuartet[3] << endl;
+                _logger->info("Searching for first of a quartet at cycle {}: {} {} {} {}",
+                              fFirstCycle, fQuartet[0], fQuartet[1], fQuartet[2], fQuartet[3]);
                 fFirstCycle++;
               }
             } else {
@@ -563,6 +539,22 @@ Int_t THcHelicity::End(THaRunBase*) {
   //  for( Int_t i = 0; i < NHIST; ++i )
   //    fHisto[i]->Write();
 
+  if (fRecommendedFreq < 0.0) {
+    fRecommendedFreq = fFreq * (1 - (fPeriodCheck - 0.5) / fCycle);
+  }
+  if (TMath::Abs(1 - fRecommendedFreq / fFreq) >= 0.5e-6) {
+    cout << "------------- HELICITY DECODING ----------------------" << endl;
+    cout << "Actual helicity reversal frequency differs from \"helicity_freq\" value" << endl;
+    cout << "If there are helicity decoding errors beyond the start of the run, " << endl;
+    streamsize ss = cout.precision();
+    cout.precision(10);
+    cout << "try replacing helicity_freq value of " << fFreq << " with " << fRecommendedFreq
+         << endl;
+    cout << "If that still gives helicity errors, try " << 0.9999999 * fRecommendedFreq << endl;
+    cout.precision(ss);
+    cout << "------------------------------------------------------" << endl;
+  }
+
   return 0;
 }
 
@@ -577,6 +569,7 @@ void THcHelicity::SetDebug(Int_t level) {
 
 //_____________________________________________________________________________
 void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t missedcycles) {
+
   //  static const char* const here = "THcHelicity::LoadHelicity";
   int quartetphase = (cyclecount - fFirstCycle) % 4;
   fnQrt            = quartetphase;
@@ -618,14 +611,14 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
       if (fReportedHelicity == kUnknown) {
         fNBits             = 0;
         fRingSeed_reported = 0;
-      } else if (fNBits == fMAXBIT) {
-        cout << "Seed Found    " << bitset<32>(fRingSeed_reported) << " at cycle " << cyclecount
-             << " with first cycle " << fFirstCycle << endl;
+      } else if (fNBits == fMAXBIT + 0) {
+        _logger->info("Seed Found {:32b} at cycle {} with first cycle {}", fRingSeed_reported,
+                      cyclecount, fFirstCycle);
         if (fglHelicityScaler) {
-          cout << "Scaler Seed " << bitset<32>(fScalerSeed) << endl;
+          _logger->info("Scaler Seed {:32b}", fScalerSeed);
         }
         Int_t backseed = GetSeed30(fRingSeed_reported);
-        cout << "Seed at cycle " << fFirstCycle << " should be " << hex << backseed << dec << endl;
+        _logger->info("Seed at cycle {} should be {:#x}", fFirstCycle, backseed);
         // Create the "actual seed"
         fRingSeed_actual = fRingSeed_reported;
         for (Int_t i = 0; i < fHelDelay / 4; i++) {
@@ -633,7 +626,24 @@ void THcHelicity::LoadHelicity(Int_t reportedhelicity, Int_t cyclecount, Int_t m
         }
         fQuartetStartHelicity          = (fRingSeed_actual & 1) ? kPlus : kMinus;
         fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
-      }
+      } else if (fglHelicityScaler && fNBits > 2) { // Try the scalers
+        if (fglHelicityScaler->IsSeedGood()) {
+          Int_t scalerseed   = fglHelicityScaler->GetReportedSeed();
+          fRingSeed_reported = RanBit30(scalerseed);
+          _logger->info(" -- Getting seed from scalers -- ");
+          cout << " Seed Found   " << bitset<32>(fRingSeed_reported) << " at cycle " << cyclecount
+               << " with first cycle " << fFirstCycle << endl;
+          cout << "Scaler Seed  " << bitset<32>(scalerseed) << endl;
+          // Create the "actual seed"
+          fRingSeed_actual = fRingSeed_reported;
+          for (Int_t i = 0; i < fHelDelay / 4; i++) {
+            fRingSeed_actual = RanBit30(fRingSeed_actual);
+          }
+          fQuartetStartHelicity          = (fRingSeed_actual & 1) ? kPlus : kMinus;
+          fQuartetStartPredictedHelicity = (fRingSeed_reported & 1) ? kPlus : kMinus;
+          fNBits                         = fMAXBIT + 0;
+        }
+      } 
       fActualHelicity = kUnknown;
     } // Need to change this to build seed even when not at start of quartet
   } else {

src/THcHelicity.h

@@ -9,9 +9,10 @@
 //
 ////////////////////////////////////////////////////////////////////////
 
+#include "hcana/Logger.h"
+
 #include "THaHelicityDet.h"
 #include "THcHelicityReader.h"
-#include "hcana/Logger.h"
 
 class TH1F;
 class THcHelicityScaler;
@@ -49,10 +50,12 @@ protected:
   Int_t    fFirstCycle;
   Bool_t   fFixFirstCycle;
   Double_t fFreq;
+  Double_t fRecommendedFreq;
 
   Double_t fTIPeriod; // Reversal period in TI time units
 
   Double_t fPeriodCheck;
+  Double_t fPeriodCheckOffset;
   Double_t fCycle;
 
   Bool_t   fFirstEvProcessed;

src/THcHelicityReader.cxx

@@ -45,9 +45,10 @@ Int_t THcHelicityReader::ReadDatabase(const char* /*dbfilename*/, const char* /*
                                       const TDatime& /*date*/, int /*debug_flag*/) {
 
   // Eventually get these from the parameter file
+  const static char* const here = "THcHelicityReader::ReadDatabase";
 
   // SHMS settings see https://logbooks.jlab.org/entry/3614445
-  cout << "THcHelicityReader: Helicity information from ROC 2 (SHMS)" << endl;
+  _logger->info("{}: Helicity information from ROC 2 (SHMS)", here);
   SetROCinfo(kHel, 2, 14, 9);
   SetROCinfo(kHelm, 2, 14, 8);
   SetROCinfo(kMPS, 2, 14, 10);
@@ -106,8 +107,7 @@ Int_t THcHelicityReader::ReadData(const THaEvData& evdata) {
   // Check if ROC info is correct
 
   if (!evdata.GetModule(fROCinfo[kTime].roc, fROCinfo[kTime].slot)) {
-    cout << "THcHelicityReader: ROC 2 not found" << endl;
-    cout << "Changing to ROC 1 (HMS)" << endl;
+    _logger->info("{}: ROC 2 not found, changing to ROC 1 (HMS)", here);
     SetROCinfo(kHel, 1, 18, 9);
     SetROCinfo(kHelm, 1, 18, 8);
     SetROCinfo(kMPS, 1, 18, 10);
@@ -122,8 +122,7 @@ Int_t THcHelicityReader::ReadData(const THaEvData& evdata) {
       (UInt_t)evdata.GetData(fROCinfo[kTime].roc, fROCinfo[kTime].slot, fROCinfo[kTime].index, 0);
   // Check again if ROC info is correct
   if (titime == 0 && fTITime_last == 0) {
-    cout << "THcHelicityReader: ROC 2 not found" << endl;
-    cout << "Changing to ROC 1 (HMS)" << endl;
+    _logger->info("{}: ROC 2 not found, changing to ROC 1 (HMS)", here);
     SetROCinfo(kHel, 1, 18, 9);
     SetROCinfo(kHelm, 1, 18, 8);
     SetROCinfo(kMPS, 1, 18, 10);

src/THcHelicityScaler.h

@@ -9,89 +9,171 @@
 /////////////////////////////////////////////////////////////////////
 
 #include "THaEvtTypeHandler.h"
+#include "THcScalerEvtHandler.h"
+
 #include "Decoder.h"
-#include <string>
-#include <vector>
-#include <algorithm>
-#include <set>
-#include "TTree.h"
+#include "THaRunBase.h"
 #include "TString.h"
+#include "TTree.h"
+#include "hcana/Logger.h"
+#include <algorithm>
 #include <cstring>
+#include <nlohmann/json.hpp>
+#include <set>
+#include <string>
+#include <vector>
 
 class THcHelicity;
+class HCScalerLoc;
 
-class THcHelicityScaler : public THaEvtTypeHandler {
+class THcHelicityScaler : public hcana::ConfigLogging<THaEvtTypeHandler> {
 
 public:
-
   THcHelicityScaler(const char*, const char*);
   virtual ~THcHelicityScaler();
 
-  Int_t Analyze(THaEvData *evdata);
-  Int_t AnalyzeBuffer(UInt_t *rdata);
-  Int_t AnalyzeHelicityScaler(UInt_t *p);
-
-  virtual EStatus Init( const TDatime& run_time);
-  virtual Int_t   ReadDatabase(const TDatime& date );
-  virtual Int_t End( THaRunBase* r=0 );
-
-  virtual void SetUseFirstEvent(Bool_t b = kFALSE) {fUseFirstEvent = b;}
-  virtual void SetDelayedType(int evtype);
-  virtual void SetROC(Int_t roc) {fROC=roc;}
-  virtual void SetBankID(Int_t bankid) {fBankID=bankid;}
-  virtual void SetNScalerChannels(Int_t n) {fNScalerChannels = n;}
-  virtual Int_t GetNevents() { return fNTrigsInBuf;}
-  virtual Int_t GetNcycles() { return fNTriggers;}
-  virtual Int_t GetEvNum() { return evNumber;}
-  virtual Int_t* GetHelicityHistoryP() {return fHelicityHistory;}
-  virtual Int_t GetReportedSeed() {return fRingSeed_reported;}
-  virtual Int_t GetReportedActual() {return fRingSeed_actual;}
-  virtual Bool_t IsSeedGood() {return fNBits>=30;}
+  Int_t           Analyze(THaEvData* evdata);
+  Int_t           AnalyzeBuffer(UInt_t* rdata);
+  Int_t           AnalyzeHelicityScaler(UInt_t* p);
+  virtual EStatus Init(const TDatime& run_time);
+  virtual Int_t   ReadDatabase(const TDatime& date);
+  virtual Int_t   End(THaRunBase* r = 0);
+
+  virtual void   SetUseFirstEvent(Bool_t b = kFALSE) { fUseFirstEvent = b; }
+  virtual void   SetDelayedType(int evtype);
+  virtual void   SetROC(Int_t roc) { fROC = roc; }
+  virtual void   SetBankID(Int_t bankid) { fBankID = bankid; }
+  virtual void   SetNScalerChannels(Int_t n) { fNScalerChannels = n; }
+  virtual Int_t  GetNevents() { return fNTrigsInBuf; }
+  virtual Int_t  GetNcycles() { return fNTriggers; }
+  virtual Int_t  GetEvNum() { return evNumber; }
+  virtual Int_t* GetHelicityHistoryP() { return fHelicityHistory; }
+  virtual Int_t  GetReportedSeed() { return fRingSeed_reported; }
+  virtual Int_t  GetReportedActual() { return fRingSeed_actual; }
+  virtual Bool_t IsSeedGood() { return fNBits >= 30; }
+
+  Double_t GetPlusCharge(const std::string& name);
+  Double_t GetMinusCharge(const std::string& name);
+
+  // utility function to write out json helicity file
+  void WriteJson(const std::string& path) const;
 
 private:
+  //------------C.Y. Sep 20, 2020 :: Added Utility Function Prototypes----------------
+  void          AddVars(TString name, TString desc, UInt_t iscal, UInt_t ichan, UInt_t ikind);
+  void          DefVars();
+  static size_t FindNoCase(const std::string& sdata, const std::string& skey);
+
+  std::vector<Decoder::GenScaler*> scalers;
+  std::vector<HCScalerLoc*>        scalerloc;
 
   Int_t RanBit30(Int_t ranseed);
-  void MakeParms();
+  void  MakeParms();
 
   UInt_t fBankID;
   // Helicity Scaler variables
-  Int_t fNTrigsInBuf;		/* # of helicity scaler reads in last event */
+  Int_t fNTrigsInBuf; /* # of helicity scaler reads in last event */
   Int_t fNTriggers;
   Int_t fFirstCycle;
   Int_t fHelicityHistory[200];
   //
   Bool_t fUseFirstEvent;
-  Int_t fDelayedType;
+  Int_t  fDelayedType;
 
   Int_t fRingSeed_reported;
   Int_t fRingSeed_actual;
   Int_t fNBits;
 
-  Int_t fNTriggersPlus;
-  Int_t fNTriggersMinus;
+  Int_t                 fNTriggersPlus;
+  Int_t                 fNTriggersMinus;
   std::vector<Double_t> fHScalers[2];
-  Int_t fGateCount[2];
+  Int_t                 fGateCount[2];
   std::vector<Double_t> fScalerSums, fAsymmetry, fAsymmetryError;
-  std::vector<Double_t> fCharge, fChargeAsymmetry;
-  Double_t fTime;
+  std::vector<Double_t> fCharge;
+  Double_t              fTime;
+  Double_t              fTimePlus;
+  Double_t              fTimeMinus;
+  Double_t              fTriggerAsymmetry;
+
+  //---- C.Y.: 12/14/2020  Variables for quartet-by-quartet asymmetry/error calculations ----
+  Bool_t fHaveCycle[4];
+
+  Int_t fQuartetCount; // keep track of number of quartets
+
+  // quartet-by-quartet time asymmetry variables
+  Double_t fTimeCycle[4];
+  Double_t fTimeSum;
   Double_t fTimeAsymmetry;
-  Double_t fTriggerAsymmetry;
+  Double_t fTimeAsymmetryError;
+  Double_t fTimeAsymSum;
+  Double_t fTimeAsymSum2;
+
+  // quartet-by-quartet scaler counts asymmetry variables
+  std::vector<Double_t> fScalCycle[4];
+  std::vector<Double_t> fScalSum; // reminder: need to initialize
+  std::vector<Double_t> fScalAsymmetry;
+  std::vector<Double_t> fScalAsymmetryError;
+  std::vector<Double_t> fScalAsymSum;
+  std::vector<Double_t> fScalAsymSum2;
+
+  // quartet-by-quartet charge asymmetry variables
+  std::vector<Double_t> fChargeCycle[4];
+  std::vector<Double_t> fChargeSum;
+  std::vector<Double_t> fChargeAsymmetry;
+  std::vector<Double_t> fChargeAsymmetryError;
+  std::vector<Double_t> fChargeAsymSum;
+  std::vector<Double_t> fChargeAsymSum2;
+
+  //----------------------
+
+  //----C.Y. Nov 26, 2020----
+  std::vector<Double_t> fScalerChan;
 
   std::vector<UInt_t*> fDelayedEvents;
-  Int_t fROC;
-  Int_t fNScalerChannels;	// Number of scaler channels/event
-
-  Int_t fNumBCMs;
-  std::vector<Double_t> fBCM_Gain, fBCM_Offset;
-  std::vector <std::string> fBCM_Name;
+  Int_t                fROC;
+  Int_t                fNScalerChannels; // Number of scaler channels/event
+
+  Int_t                    fNumBCMs;
+  std::vector<Double_t>    fBCM_Gain, fBCM_Offset;
+  std::vector<std::string> fBCM_Name;
+  //---C.Y. Sep 2020 : Added additional BCM-related variables--
+  std::vector<Double_t> fBCM_SatOffset;
+  std::vector<Double_t> fBCM_SatQuadratic;
+  std::vector<Double_t> fBCM_delta_charge;
+  Double_t              fTotalTime;
+  Double_t              fDeltaTime;
+  Double_t              fPrevTotalTime;
+  Double_t              fbcm_Current_Threshold;
+  Double_t              fClockFreq;
+  Int_t                 fbcm_Current_Threshold_Index;
+
+  //----C.Y. Sep 20, 2020 : Added additional variables-----
+  // (required by utility functions and scaler tree output)
+  UInt_t                evcount;
+  Double_t              evcountR;
+  UInt_t                evNumber;
+  Double_t              evNumberR;
+  Double_t              actualHelicityR;
+  Double_t              quartetphaseR;
+  Int_t                 Nvars, ifound, fNormIdx, fNormSlot, nscalers;
+  std::vector<Double_t> dvars;
+  std::vector<Double_t> dvarsFirst;
+  TTree*                fScalerTree;
+  Bool_t                fOnlySyncEvents;
+  Bool_t                fOnlyBanks;
+  Int_t                 fClockChan;
+  UInt_t                fLastClock;
+  std::set<UInt_t>      fRocSet;
+  std::set<UInt_t>      fModuleSet;
+  //--------------------------------------------------------
+
+  // json file with helicity info
+  nlohmann::json j_helicity;
 
   THcHelicityScaler(const THcHelicityScaler& fh);
   THcHelicityScaler& operator=(const THcHelicityScaler& fh);
 
-  UInt_t evNumber;
-
-  ClassDef(THcHelicityScaler,0)  // Scaler Event handler
-
+  ClassDef(THcHelicityScaler, 0) // Scaler Event handler
 };
 
 #endif

src/THcHitList.cxx

@@ -263,10 +263,14 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	fRefIndexMaps[i].hashit = kFALSE;
 	Bool_t goodreftime=kFALSE;
 	Int_t reftime = 0;
+	Int_t prevtime = 0;
+	Int_t difftime = 0;
 	for(Int_t ihit=0; ihit<nrefhits; ihit++) {
 	  reftime = evdata.GetData(Decoder::kPulseTime,fRefIndexMaps[i].crate,
 				   fRefIndexMaps[i].slot, fRefIndexMaps[i].channel,ihit);
 	  reftime += 64*timeshift;
+	  if (ihit != 0) difftime=reftime-prevtime;
+	  prevtime = reftime;
 	  if(reftime >= fADC_RefTimeCut) {
 	    goodreftime = kTRUE;
 	    break;
@@ -274,6 +278,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	}
 	if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
 	  fRefIndexMaps[i].reftime = reftime;
+	  fRefIndexMaps[i].refdifftime = difftime;
 	  fRefIndexMaps[i].hashit = kTRUE;
 	}
       } else {			// Assume this is a TDC
@@ -285,9 +290,13 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	// then fTDC_RefTimeCut
 	Bool_t goodreftime=kFALSE;
 	Int_t reftime = 0;
+	Int_t prevtime = 0;
+	Int_t difftime = 0;
 	for(Int_t ihit=0; ihit<nrefhits; ihit++) {
 	  reftime = evdata.GetData(fRefIndexMaps[i].crate,fRefIndexMaps[i].slot,
 				   fRefIndexMaps[i].channel,ihit);
+	  if( ihit != 0) difftime=reftime-prevtime;
+	    prevtime=reftime;
 	  if(reftime >= fTDC_RefTimeCut) {
 	    goodreftime = kTRUE;
 	    break;
@@ -295,6 +304,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	}
 	if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
 	    fRefIndexMaps[i].reftime = reftime;
+	    fRefIndexMaps[i].refdifftime = difftime;
 	    fRefIndexMaps[i].hashit = kTRUE;
 	}
       }
@@ -360,8 +370,12 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  Int_t nrefhits = evdata.GetNumHits(d->crate,d->slot,d->refchan);
 	  Bool_t goodreftime=kFALSE;
 	  Int_t reftime=0;
+	  Int_t prevtime=0;
+	  Int_t difftime=0;
 	  for(Int_t ihit=0; ihit<nrefhits; ihit++) {
 	    reftime = evdata.GetData(d->crate, d->slot, d->refchan, ihit);
+	    if (ihit != 0 ) difftime=reftime-prevtime;
+              prevtime = reftime;
 	    if(reftime >= fTDC_RefTimeCut) {
 	      goodreftime = kTRUE;
 	      break;
@@ -371,6 +385,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  // hits make the RefTimeCut
 	  if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
 	    rawhit->SetReference(signal, reftime);
+	    rawhit->SetReferenceDiff(signal, difftime);
 	  } else if (!suppresswarnings) {
 	    cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
 	      " missing for (" << d->crate << ", " << d->slot <<
@@ -381,6 +396,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  if(d->refindex >=0 && d->refindex < fNRefIndex) {
 	    if(fRefIndexMaps[d->refindex].hashit) {
 	      rawhit->SetReference(signal, fRefIndexMaps[d->refindex].reftime);
+	      rawhit->SetReferenceDiff(signal, fRefIndexMaps[d->refindex].refdifftime);
 	    } else {
 	      if(!suppresswarnings) {
 		cout << "HitList(event=" << evdata.GetEvNum() << "): refindex " << d->refindex <<
@@ -451,6 +467,8 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 					       d->crate, d->slot, d->refchan);
 	  Bool_t goodreftime=kFALSE;
 	  Int_t reftime = 0;
+	  Int_t prevtime = 0;
+	  Int_t difftime = 0;
 	  timeshift=0;
 	  if(fTISlot>0) {		// Get the trigger time for this module
 	    if(fTrigTimeShiftMap.find(d->slot)
@@ -465,6 +483,8 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  for(Int_t ihit=0; ihit<nrefhits; ihit++) {
 	    reftime = evdata.GetData(Decoder::kPulseTime, d->crate, d->slot, d->refchan, ihit);
 	    reftime += 64*timeshift;
+	    if (ihit != 0) difftime=reftime-prevtime;
+	    prevtime=reftime;
 	    if(reftime >= fADC_RefTimeCut) {
 	      goodreftime=kTRUE;
 	      break;
@@ -474,6 +494,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  // hits make the RefTimeCut
 	  if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
 	    rawhit->SetReference(signal, reftime);
+	    rawhit->SetReferenceDiff(signal, difftime);
 	  } else if (!suppresswarnings) {
 #ifndef SUPPRESSMISSINGADCREFTIMEMESSAGES
 	    cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
@@ -486,6 +507,7 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  if(d->refindex >=0 && d->refindex < fNRefIndex) {
 	    if(fRefIndexMaps[d->refindex].hashit) {
 	      rawhit->SetReference(signal, fRefIndexMaps[d->refindex].reftime);
+	      rawhit->SetReferenceDiff(signal, fRefIndexMaps[d->refindex].refdifftime);
 	    } else {
 	      if(!suppresswarnings) {
 #ifndef SUPPRESSMISSINGADCREFTIMEMESSAGES

src/THcHitList.h

@@ -73,6 +73,7 @@ protected:
     Int_t slot;
     Int_t channel;
     Int_t reftime;
+    Int_t refdifftime;
   };
   std::vector<RefIndexMap> fRefIndexMaps;
   // Should this be a sparse list instead in case user