From cd9cc78cc34a463907fee7c0103aef1b4092beb0 Mon Sep 17 00:00:00 2001
From: Chao Peng <thor1009@gmail.com>
Date: Tue, 31 Dec 2019 20:58:35 -0600
Subject: [PATCH] Merge Steve's Helicity Scaler code into the repo

---
 src/THcHelicity.cxx       |   2 +
 src/THcHelicityScaler.cxx | 522 +++++++++++++++++++++++++++-----------
 src/THcHelicityScaler.h   | 104 ++++----
 3 files changed, 441 insertions(+), 187 deletions(-)

diff --git a/src/THcHelicity.cxx b/src/THcHelicity.cxx
index 0b36432..a242703 100644
--- a/src/THcHelicity.cxx
+++ b/src/THcHelicity.cxx
@@ -270,6 +270,7 @@ 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();
@@ -303,6 +304,7 @@ Int_t THcHelicity::Decode(const THaEvData& evdata) {
       }
     }
   }
+#endif
 
   if (fHelDelay == 0) { // If no delay actual=reported (but zero if in MPS)
     fActualHelicity = fIsMPS ? kUnknown : fReportedHelicity;
diff --git a/src/THcHelicityScaler.cxx b/src/THcHelicityScaler.cxx
index 8f4ce39..2ce79bd 100644
--- a/src/THcHelicityScaler.cxx
+++ b/src/THcHelicityScaler.cxx
@@ -3,96 +3,193 @@
 
 \brief Event handler for Hall C helicity scalers
 
-Scalers not yet implemented.  For now just picks the helicity control
-bits out of the scaler words
-
-~~~
-     gHaEvtHandlers->Add (new THcHelicityScaler("H","HC helicity scalers"));
 ~~~
-To enable debugging you may try this in the setup script
 ~~~
      THcHelcityScaler *hhelscaler = new THcHelicityScaler("H","HC helicity scalers");
-     hscaler->SetDebugFile("HHelScaler.txt");
+     // hscaler->SetDebugFile("HHelScaler.txt");
+     hhelscaler->SetROC(8);   // 5 for HMS defaults to 8 for SHMS
+     hhelscaler->SetBankID(0x9801); // Will default to this
      gHaEvtHandlers->Add (hhelscaler);
 ~~~
 \author
 */
-#include <cstdio>
-#include <cstdlib>
-#include <cstring>
-#include <iostream>
-#include <iterator>
-#include <map>
-#include <sstream>
 
-#include "Helper.h"
+//#include "THaEvtTypeHandler.h"
+#include "THcHelicityScaler.h"
 #include "THaCodaData.h"
 #include "THaEvData.h"
-#include "THaEvtTypeHandler.h"
-#include "THaGlobals.h"
-#include "THaVarList.h"
 #include "THcGlobals.h"
-#include "THcHelicity.h"
-#include "THcHelicityScaler.h"
+#include "THaGlobals.h"
 #include "THcParmList.h"
-#include "TMath.h"
+#include "THcHelicity.h"
 #include "TNamed.h"
-#include "TROOT.h"
+#include "TMath.h"
 #include "TString.h"
+#include "TROOT.h"
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+#include <sstream>
+#include <map>
+#include <bitset>
+#include <iterator>
+#include "THaVarList.h"
 #include "VarDef.h"
+#include "Helper.h"
 
 using namespace std;
 using namespace Decoder;
 
 static const UInt_t ICOUNT    = 1;
 static const UInt_t IRATE     = 2;
-static const UInt_t ICURRENT  = 3;
+static const UInt_t ICURRENT = 3;
 static const UInt_t ICHARGE   = 4;
-static const UInt_t ITIME     = 5;
-static const UInt_t ICUT      = 6;
+static const UInt_t ITIME   = 5;
+static const UInt_t ICUT = 6;
 static const UInt_t MAXCHAN   = 32;
 static const UInt_t defaultDT = 4;
 
-THcHelicityScaler::THcHelicityScaler(const char* name, const char* description)
-    : THaEvtTypeHandler(name, description), fBCM_Gain(0), fBCM_Offset(0), fBCM_delta_charge(0),
-      fBankID(9801), evcount(0), evcountR(0.0), ifound(0), fUseFirstEvent(kTRUE),
-      fOnlySyncEvents(kFALSE), fOnlyBanks(kFALSE), fDelayedType(-1), fClockChan(-1), fLastClock(0),
-      fClockOverflows(0) {
-  fRocSet.clear();
+THcHelicityScaler::THcHelicityScaler(const char *name, const char* description)
+  : THaEvtTypeHandler(name,description),
+    fBankID(9801),
+    fUseFirstEvent(kTRUE),
+    fDelayedType(-1),
+    fBCM_Gain(0), fBCM_Offset(0)
+{
+  fROC=-1;
+  fNScalerChannels = 32;
+
+  AddEvtType(1);
+  AddEvtType(2);
+  AddEvtType(4);
+  AddEvtType(5);
+  AddEvtType(6);
+  AddEvtType(7);
+  SetDelayedType(129);
+
 }
 
-THcHelicityScaler::~THcHelicityScaler() {
-  delete[] fBCM_Gain;
-  delete[] fBCM_Offset;
-  delete[] fBCM_delta_charge;
+THcHelicityScaler::~THcHelicityScaler()
+{
+  delete [] fBCM_Gain;
+  delete [] fBCM_Offset;
 
-  for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
-    delete[] * it;
+  for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
+       it != fDelayedEvents.end(); ++it )
+    delete [] *it;
   fDelayedEvents.clear();
 }
 
-Int_t THcHelicityScaler::End(THaRunBase*) {
+Int_t THcHelicityScaler::End( THaRunBase* )
+{
   // Process any delayed events in order received
 
-  cout << "THcHelicityScaler::End Analyzing " << fDelayedEvents.size() << " delayed scaler events"
-       << endl;
-  for (std::vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end();
-       ++it) {
+  cout << "THcHelicityScaler::End Analyzing " << fDelayedEvents.size() << " delayed helicity scaler events" << endl;
+    for(std::vector<UInt_t*>::iterator it = fDelayedEvents.begin();
+      it != fDelayedEvents.end(); ++it) {
     UInt_t* rdata = *it;
-    AnalyzeBuffer(rdata, kFALSE);
+    AnalyzeBuffer(rdata);
   }
 
-  for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
-    delete[] * it;
+  for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
+       it != fDelayedEvents.end(); ++it )
+    delete [] *it;
   fDelayedEvents.clear();
 
+  //  cout << " -- Helicity Scalers -- " << endl;
+  for(Int_t i=0;i<fNScalerChannels;i++) {
+    if(fScalerSums[i]>0.5) {
+      fAsymmetry[i] = (fHScalers[0][i]-fHScalers[1][i])/fScalerSums[i];
+      fAsymmetryError[i] = 2*TMath::Sqrt(fHScalers[0][i]*fHScalers[1][i]
+                                        *fScalerSums[i])
+        /(fScalerSums[i]*fScalerSums[i]);
+    } else {
+      fAsymmetry[i] = 0.0;
+      fAsymmetryError[i] = 0.0;
+    }
+    //    printf("%2d %12.0f %12.0f %12.0f %12.8f\n",i,fScalerSums[i],
+    //           fHScalers[0][i],fHScalers[1][i],
+    //           fAsymmetry[i]);
+  }
+  //  cout << " ---------------------- " << endl;
+
+  // Compute Charge Asymmetries
+  std::map<std::string, Int_t> bcmindex;
+  bcmindex["BCM1"] = 0;
+  bcmindex["BCM2"] = 2;
+  //  bcmindex["Unser"] = 6;
+  bcmindex["BCM4A"] = 10;
+  bcmindex["BCM4B"] = 4;
+  bcmindex["BCM4C"] = 12;
+  //  bcmindex["1MHz"] = 8;
+  Int_t clockindex=8;
+  Double_t clockfreq=1000000.0;
+  Double_t pclock = fHScalers[0][clockindex];
+  Double_t mclock = fHScalers[1][clockindex];
+  cout << " -- Beam Charge Asymmetries -- " << endl;
+  for(Int_t i=0;i<fNumBCMs;i++) {
+    if(bcmindex.find(fBCM_Name[i]) != bcmindex.end()) {
+      Int_t index=bcmindex[fBCM_Name[i]];
+      Double_t pcounts = fHScalers[0][index];
+      Double_t mcounts = fHScalers[1][index];
+      //      cout << index << " " << fBCM_Name[i] << " " << pcounts << " " << mcounts
+      //           << " " << fBCM_Gain[i]
+      //                 << " " << fBCM_Offset[i] << endl;
+      Double_t pcharge = (pcounts - (pclock/clockfreq)*fBCM_Offset[i])
+        /fBCM_Gain[i];
+      Double_t mcharge = (mcounts - (mclock/clockfreq)*fBCM_Offset[i])
+        /fBCM_Gain[i];
+      fCharge[i] = pcharge+mcharge;
+      if(fCharge[i]>0.0) {
+        fChargeAsymmetry[i] = (pcharge-mcharge)/fCharge[i];
+      } else {
+        fChargeAsymmetry[i] = 0.0;
+      }
+      printf("%6s %12.2f %12.8f\n",fBCM_Name[i].c_str(),fCharge[i],fChargeAsymmetry[i]);
+    }
+  }
+  fTime = (pclock+mclock)/clockfreq;
+  if(pclock+mclock>0) {
+    fTimeAsymmetry = (pclock-mclock)/(pclock+mclock);
+  } else {
+    fTimeAsymmetry = 0.0;
+  }
+  printf("TIME(s)%12.2f %12.8f\n",fTime,fTimeAsymmetry);
+  if(fNTriggersPlus+fNTriggersMinus > 0) {
+    fTriggerAsymmetry = ((Double_t) (fNTriggersPlus-fNTriggersMinus))/(fNTriggersPlus+fNTriggersMinus);
+  } else {
+    fTriggerAsymmetry = 0.0;
+  }
+  cout << " ----------------------------- " << endl;
   return 0;
 }
 
-Int_t THcHelicityScaler::ReadDatabase(const TDatime& date) {
+
+Int_t THcHelicityScaler::ReadDatabase(const TDatime& date )
+{
   char prefix[2];
-  prefix[0] = 'g';
-  prefix[1] = '\0';
+  prefix[0]='g'; prefix[1]='\0';
+
+  fNumBCMs = 0;
+  string bcm_namelist;
+  DBRequest list[]={
+                    {"NumBCMs",&fNumBCMs, kInt, 0, 1},
+                    {"BCM_Names",     &bcm_namelist,       kString},
+                    {0}
+  };
+  gHcParms->LoadParmValues((DBRequest*)&list, prefix);
+  if(fNumBCMs > 0) {
+    fBCM_Gain = new Double_t[fNumBCMs];
+    fBCM_Offset = new Double_t[fNumBCMs];
+    DBRequest list2[]={
+      {"BCM_Gain",      fBCM_Gain,         kDouble, (UInt_t) fNumBCMs},
+      {"BCM_Offset",     fBCM_Offset,       kDouble,(UInt_t) fNumBCMs},
+      {0}
+    };
+    gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
+    fBCM_Name = vsplit(bcm_namelist);
+  }
 
   return kOK;
 }
@@ -108,70 +205,66 @@ void THcHelicityScaler::SetDelayedType(int evtype) {
   fDelayedType = evtype;
 }
 
-Int_t THcHelicityScaler::Analyze(THaEvData* evdata) {
+Int_t THcHelicityScaler::Analyze(THaEvData *evdata)
+{
 
-  if (!IsMyEvent(evdata->GetEvType()))
-    return -1;
+  if ( !IsMyEvent(evdata->GetEvType()) ) return -1;
 
   if (fDebugFile) {
-    *fDebugFile << endl << "---------------------------------- " << endl << endl;
-    *fDebugFile << "\nEnter THcHelicityScaler  for fName = " << fName << endl;
+    *fDebugFile << endl << "---------------------------------- "<<endl<<endl;
+    *fDebugFile << "\nEnter THcHelicityScaler  for fName = "<<fName<<endl;
     EvDump(evdata);
   }
 
-  UInt_t* rdata = (UInt_t*)evdata->GetRawDataBuffer();
+  UInt_t *rdata = (UInt_t*) evdata->GetRawDataBuffer();
 
-  if (evdata->GetEvType() == fDelayedType) { // Save this event for processing later
-    Int_t   evlen    = evdata->GetEvLength();
-    UInt_t* datacopy = new UInt_t[evlen];
+  if(evdata->GetEvType() == fDelayedType) { // Save this event for processing later
+    Int_t evlen = evdata->GetEvLength();
+    UInt_t *datacopy = new UInt_t[evlen];
     fDelayedEvents.push_back(datacopy);
-    memcpy(datacopy, rdata, evlen * sizeof(UInt_t));
+    memcpy(datacopy,rdata,evlen*sizeof(UInt_t));
     return 1;
-  } else { // A normal event
-    if (fDebugFile)
-      *fDebugFile << "\n\nTHcHelicityScaler :: Debugging event type " << dec << evdata->GetEvType()
-                  << " event num = " << evdata->GetEvNum() << endl
-                  << endl;
-    evNumber = evdata->GetEvNum();
+  } else {                         // A normal event
+    if (fDebugFile) *fDebugFile<<"\n\nTHcHelicityScaler :: Debugging event type "<<dec<<evdata->GetEvType()<< " event num = " << evdata->GetEvNum() << endl<<endl;
+    evNumber=evdata->GetEvNum();
     Int_t ret;
-    if ((ret = AnalyzeBuffer(rdata, fOnlySyncEvents))) {
+    if((ret=AnalyzeBuffer(rdata))) {
       //
     }
     return ret;
   }
+
 }
-Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync) {
-  fNevents = 0;
+Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata)
+{
+  fNTrigsInBuf = 0;
 
   // Parse the data, load local data arrays.
-  UInt_t* p = (UInt_t*)rdata;
+  UInt_t *p = (UInt_t*) rdata;
 
-  UInt_t* plast = p + *p; // Index to last word in the bank
-  Int_t   roc   = -1;
-  Int_t   evlen = *p + 1;
+  UInt_t *plast = p+*p;                // Index to last word in the bank
+  Int_t roc = -1;
+  Int_t evlen = *p+1;
 
-  ifound = 0;
+  Int_t ifound=0;
 
-  while (p < plast) {
+  while(p<plast) {
     Int_t banklen = *p;
-    p++; // point to header
+    p++;                          // point to header
 
     if (fDebugFile) {
-      *fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p - 1) << endl;
+      *fDebugFile << "Bank: " << hex << *p << dec << " len: " << *(p-1) << endl;
     }
-    if ((*p & 0xff00) == 0x1000) { // Bank Containing banks
-      if (evlen - *(p - 1) > 1) {  // Don't use overall event header
-        roc = (*p >> 16) & 0xf;
-        if (fDebugFile)
-          *fDebugFile << "ROC: " << roc << " " << evlen << " " << *(p - 1) << hex << " " << *p
-                      << dec << endl;
-        //		cout << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p <<
-        //dec << endl;
-        if (fRocSet.find(roc) == fRocSet.end()) { // Not a ROC with helicity scaler
-          p += *(p - 1) - 1;                      // Skip to next ROC
+    if((*p & 0xff00) == 0x1000) {        // Bank Containing banks
+      if(evlen-*(p-1) > 1) { // Don't use overall event header
+        roc = (*p>>16) & 0xf;
+        if(fDebugFile) *fDebugFile << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p << dec << endl;
+//                cout << "ROC: " << roc << " " << evlen << " " << *(p-1) << hex << " " << *p << dec << endl;
+        if(roc != fROC) { // Not a ROC with helicity scaler
+          p+=*(p-1)-1;                // Skip to next ROC
         }
       }
-      p++; // Now pointing to a bank in the bank
+      p++;                                // Now pointing to a bank in the bank
     } else if (((*p & 0xff00) == 0x100) && (*p != 0xC0000100)) {
       // Bank containing integers.  Look for scalers
       // This is either ROC bank containing integers or
@@ -180,115 +273,258 @@ Int_t THcHelicityScaler::AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync) {
       // Assume that very first word is a scaler header
       // At any point in the bank where the word is not a matching
       // header, we stop.
-      UInt_t tag = (*p >> 16) & 0xffff; // Bank ID (ROC #)
-                                        //          UInt_t num = (*p) & 0xff;
-      UInt_t* pnext = p + banklen;      // Next bank
-      p++;                              // First data word
+      UInt_t tag = (*p>>16) & 0xffff; // Bank ID (ROC #)
+  //          UInt_t num = (*p) & 0xff;
+      UInt_t *pnext = p+banklen;        // Next bank
+      p++;                        // First data word
       // If the bank is not a helicity scaler bank
       // or it is not one of the ROC containing helcity scaler data
       // skip to the next bank
-      // cout << "BankID=" << tag << endl;
+      //cout << "BankID=" << tag << endl;
 
       if (tag != fBankID) {
-        p = pnext; // Fall through to end of the above else if
-                   //	cout << "  Skipping to next bank" << endl;
+        p = pnext;                // Fall through to end of the above else if
+        //        cout << "  Skipping to next bank" << endl;
 
       } else {
         // This is a helicity scaler bank
-        if (roc == 5) {
-          Int_t nevents = (banklen - 2) / 32;
-          // cout << "# of helicity events in bank:" << " " << nevents << endl;
+        if (roc == fROC) {
+          Int_t nevents = (banklen-2)/fNScalerChannels;
+          //cout << "# of helicity events in bank:" << " " << nevents << endl;
           if (nevents > 100) {
             cout << "Error! Beam off for too long" << endl;
           }
-          fNevents = nevents;
-          fNcycles += nevents;
 
+          fNTrigsInBuf = 0;
           // Save helcitiy and quad info for THcHelicity
-          for (Int_t iev = 0; iev < nevents; iev++) { // find number of helicity events in each bank
-            Int_t nentries        = 32 * iev + 2;
-            fHelicityHistory[iev] = (p[nentries - 1] >> 30) & 0x3;
-            //	    cout << "H: " << evNumber << endl;
+          for (Int_t iev = 0; iev < nevents; iev++) {  // find number of helicity events in each bank
+            Int_t index = fNScalerChannels*iev+1;
+            AnalyzeHelicityScaler(p+index);
+            //            cout << "H: " << evNumber << endl;
           }
         }
       }
 
-      while (p < pnext) {
+      while(p < pnext) {
         Int_t nskip = 0;
-        if (fDebugFile) {
+        if(fDebugFile) {
           *fDebugFile << "Scaler Header: " << hex << *p << dec;
         }
-        if (nskip == 0) {
-          if (fDebugFile) {
+        if(nskip == 0) {
+          if(fDebugFile) {
             *fDebugFile << endl;
           }
-          break; // Didn't find a matching header
+          break;        // Didn't find a matching header
         }
         p = p + nskip;
       }
       p = pnext;
     } else {
-      p = p + *(p - 1); // Skip to next bank
+      p = p+*(p-1);                // Skip to next bank
     }
+
   }
 
   if (fDebugFile) {
-    *fDebugFile << "Finished with decoding.  " << endl;
-    *fDebugFile << "   Found flag   =  " << ifound << endl;
+    *fDebugFile << "Finished with decoding.  "<<endl;
+    *fDebugFile << "   Found flag   =  "<<ifound<<endl;
   }
 
-  // HMS has headers which are different from SOS, but both are
-  // event type 0 and come here.  If you found no headers, return.
-
-  if (!ifound)
-    return 0;
+  if (!ifound) return 0;
 
   return 1;
+
+
 }
 
-// Int_t THcHelicityScaler::AnalyzeHelicityScaler(UInt_t *p)
-//{
-//}
+Int_t THcHelicityScaler::AnalyzeHelicityScaler(UInt_t *p)
+{
+  Int_t hbits = (p[0]>>30) & 0x3; // quartet and helcity bits in scaler word
+  Bool_t isquartet = (hbits&2) != 0;
+  Int_t ispos = hbits&1;
+  Int_t actualhelicity = 0;
+  fHelicityHistory[fNTrigsInBuf] = hbits;
+  fNTrigsInBuf++;
+  fNTriggers++;
+
+  Int_t quartetphase = (fNTriggers-fFirstCycle)%4;
+  if(fFirstCycle >= -10) {
+    if(quartetphase == 0) {
+      Int_t predicted = RanBit30(fRingSeed_reported);
+      fRingSeed_reported = ((fRingSeed_reported<<1) | ispos) & 0x3FFFFFFF;
+      // Check if ringseed_predicted agrees with reported if(fNBits>=30)
+      if(fNBits >= 30 && predicted != fRingSeed_reported) {
+        cout << "THcHelicityScaler: Helicity Prediction Failed" << endl;
+        cout << "Reported  " << bitset<32>(fRingSeed_reported) << endl;
+        cout << "Predicted " << bitset<32>(predicted) << endl;
+      }
+      fNBits++;
+      if(fNBits==30) {
+        cout << "THcHelicityScaler: A " << bitset<32>(fRingSeed_reported) <<
+          " found at cycle " << fNTriggers << endl;
+      }
+    } else if (quartetphase == 3) {
+      if(!isquartet) {
+        cout << "THcHelicityScaler: Quartet bit expected but not set (" <<
+          fNTriggers << ")" << endl;
+        fNBits = 0;
+        fRingSeed_reported = 0;
+        fRingSeed_actual = 0;
+        fFirstCycle = -100;
+      }
+    }
+  } else {                         // First cycle not yet identified
+    if(isquartet) { // Helicity and quartet signal for next set of scalers
+      fFirstCycle = fNTriggers - 3;
+      quartetphase = (fNTriggers-fFirstCycle)%4;
+      //// Helicity at start of quartet is same as last of quartet, so we can start filling the seed
+      fRingSeed_reported = ((fRingSeed_reported<<1) | ispos) & 0x3FFFFFFF;
+      fNBits++;
+      if(fNBits==30) {
+        cout << "THcHelicityScaler: B " << bitset<32>(fRingSeed_reported) <<
+          " found at cycle " << fNTriggers << endl;
+      }
+    }
+  }
+
+  if(fNBits>=30) {
+    fRingSeed_actual = RanBit30(fRingSeed_reported);
+    fRingSeed_actual = RanBit30(fRingSeed_actual);
 
-THaAnalysisObject::EStatus THcHelicityScaler::Init(const TDatime& date) {
+#define DELAY9
+#ifdef DELAY9
+    if(quartetphase == 3) {
+      fRingSeed_actual = RanBit30(fRingSeed_actual);
+      actualhelicity = (fRingSeed_actual&1)?+1:-1;
+    } else {
+      actualhelicity = (fRingSeed_actual&1)?+1:-1;
+      if(quartetphase == 0 || quartetphase == 1) {
+        actualhelicity = -actualhelicity;
+      }
+    }
+#else
+    actualhelicity = (fRingSeed_actual&1)?+1:-1;
+    if(quartetphase == 1 || quartetphase == 2) {
+      actualhelicity = -actualhelicity;
+    }
+#endif
+  } else {
+    fRingSeed_actual = 0;
+  }
+
+  if(actualhelicity!=0) {
+    Int_t hindex = (actualhelicity>0)?0:1;
+    (actualhelicity>0)?(fNTriggersPlus++):(fNTriggersMinus++);
+    for(Int_t i=0;i<fNScalerChannels;i++) {
+      Int_t count = p[i]&0xFFFFFF; // Bottom 24 bits
+      fHScalers[hindex][i] += count;
+      fScalerSums[i] += count;
+    }
+  }
+  return(0);
+}
+//_____________________________________________________________________________
+Int_t  THcHelicityScaler::RanBit30(Int_t ranseed)
+{
+
+  UInt_t bit7    = (ranseed & 0x00000040) != 0;
+  UInt_t bit28   = (ranseed & 0x08000000) != 0;
+  UInt_t bit29   = (ranseed & 0x10000000) != 0;
+  UInt_t bit30   = (ranseed & 0x20000000) != 0;
+
+  UInt_t newbit = (bit30 ^ bit29 ^ bit28 ^ bit7) & 0x1;
+
+  ranseed =  ( (ranseed<<1) | newbit ) & 0x3FFFFFFF;
+
+  return ranseed;
+
+}
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcHelicityScaler::Init(const TDatime& date)
+{
 
   ReadDatabase(date);
 
   fStatus = kOK;
 
-  for (vector<UInt_t*>::iterator it = fDelayedEvents.begin(); it != fDelayedEvents.end(); ++it)
-    delete[] * it;
+  for( vector<UInt_t*>::iterator it = fDelayedEvents.begin();
+       it != fDelayedEvents.end(); ++it )
+    delete [] *it;
   fDelayedEvents.clear();
 
-  cout << "Howdy !  We are initializing THcHelicityScaler !!   name =   " << fName << endl;
+  cout << "Howdy !  We are initializing THcHelicityScaler !!   name =   "
+        << fName << endl;
 
-  if (eventtypes.size() == 0) {
-    eventtypes.push_back(0); // Default Event Type
+  if(eventtypes.size()==0) {
+    eventtypes.push_back(0);  // Default Event Type
   }
 
-  fRocSet.insert(5); // List ROCs that have helicity scalers
-  fRocSet.insert(8); // Should make configurable
+  if(fROC < 0) {
+    fROC = 8;                        // Default to SHMS crate
+  }
+
+  fNTriggers = 0;
+  fNTrigsInBuf = 0;
+  fFirstCycle = -100;
+  fRingSeed_reported = 0;
+  fRingSeed_actual = 0;
+  fNBits = 0;
+  fNTriggersPlus = fNTriggersMinus = 0;
+  fHScalers[0] = new Double_t[fNScalerChannels];
+  fHScalers[1] = new Double_t[fNScalerChannels];
+  fScalerSums = new Double_t[fNScalerChannels];
+  fAsymmetry = new Double_t[fNScalerChannels];
+  fAsymmetryError = new Double_t[fNScalerChannels];
+  for(Int_t i=0;i<fNScalerChannels;i++) {
+    fHScalers[0][i] = 0.0;
+    fHScalers[1][i] = 0.0;
+    fScalerSums[0] = 0.0;
+    fAsymmetry[0] = 0.0;
+    fAsymmetryError[0] = 0.0;
+  }
+
+  fCharge = new Double_t[fNumBCMs];
+  fChargeAsymmetry = new Double_t[fNumBCMs];
+
+  fTime = fTimeAsymmetry = 0;
+  fTriggerAsymmetry = 0.0;
 
-  fNcycles = 0;
-  fNevents = 0;
+  MakeParms();
 
   return kOK;
 }
 
-size_t THcHelicityScaler::FindNoCase(const string& sdata, const string& skey) {
-  // Find iterator of word "sdata" where "skey" starts.  Case insensitive.
-  string sdatalc, skeylc;
-  sdatalc = "";
-  skeylc  = "";
-  for (string::const_iterator p = sdata.begin(); p != sdata.end(); ++p) {
-    sdatalc += tolower(*p);
-  }
-  for (string::const_iterator p = skey.begin(); p != skey.end(); ++p) {
-    skeylc += tolower(*p);
-  }
-  if (sdatalc.find(skeylc, 0) == string::npos)
-    return -1;
-  return sdatalc.find(skeylc, 0);
-};
+void THcHelicityScaler::MakeParms()
+{
+  /**
+     Put Various helicity scaler results in gHcParms so they can be included in results.
+  */
+  gHcParms->Define(Form("g%s_hscaler_plus[%d]",fName.Data(),fNScalerChannels),
+                   "Plus Helcity Scalers",*(fHScalers[0]));
+  gHcParms->Define(Form("g%s_hscaler_minus[%d]",fName.Data(),fNScalerChannels),
+                   "Minus Helcity Scalers",*(fHScalers[1]));
+  gHcParms->Define(Form("g%s_hscaler_sum[%d]",fName.Data(),fNScalerChannels),
+                   "Helcity Scalers Sum",*fScalerSums);
+  gHcParms->Define(Form("g%s_hscaler_asy[%d]",fName.Data(),fNScalerChannels),
+                   "Helicity Scaler Asymmetry[%d]",*fAsymmetry);
+  gHcParms->Define(Form("g%s_hscaler_asyerr[%d]",fName.Data(),fNScalerChannels),
+                   "Helicity Scaler Asymmetry Error[%d]",*fAsymmetryError);
+  gHcParms->Define(Form("g%s_hscaler_triggers",fName.Data()),
+                   "Total Helicity Scaler Triggers",fNTriggers);
+  gHcParms->Define(Form("g%s_hscaler_triggers_plus",fName.Data()),
+                   "Positive Helicity Scaler Triggers",fNTriggersPlus);
+  gHcParms->Define(Form("g%s_hscaler_triggers_minus",fName.Data()),
+                   "Negative Helicity Scaler Triggers",fNTriggersMinus);
+  gHcParms->Define(Form("g%s_hscaler_charge[%d]",fName.Data(),fNumBCMs),
+                   "Helicity Gated Charge",*fCharge);
+  gHcParms->Define(Form("g%s_hscaler_charge_asy[%d]",fName.Data(),fNumBCMs),
+                   "Helicity Gated Charge Asymmetry",*fChargeAsymmetry);
+  gHcParms->Define(Form("g%s_hscaler_time",fName.Data()),
+                   "Helicity Gated Time (sec)",fTime);
+  gHcParms->Define(Form("g%s_hscaler_time_asy",fName.Data()),
+                   "Helicity Gated Time Asymmetry",fTimeAsymmetry);
+  gHcParms->Define(Form("g%s_hscaler_trigger_asy",fName.Data()),
+                   "Helicity Trigger Asymmetry",fTriggerAsymmetry);
+}
 
 ClassImp(THcHelicityScaler)
diff --git a/src/THcHelicityScaler.h b/src/THcHelicityScaler.h
index e065257..fff1b69 100644
--- a/src/THcHelicityScaler.h
+++ b/src/THcHelicityScaler.h
@@ -8,78 +8,94 @@
 //
 /////////////////////////////////////////////////////////////////////
 
-#include "Decoder.h"
 #include "THaEvtTypeHandler.h"
-#include "TString.h"
-#include "TTree.h"
-#include <algorithm>
-#include <cstring>
-#include <set>
+#include "Decoder.h"
 #include <string>
 #include <vector>
+#include <algorithm>
+#include <set>
+#include "TTree.h"
+#include "TString.h"
+#include <cstring>
 
 class THcHelicity;
 
 class THcHelicityScaler : public THaEvtTypeHandler {
 
 public:
+
   THcHelicityScaler(const char*, const char*);
   virtual ~THcHelicityScaler();
 
-  Int_t Analyze(THaEvData* evdata);
-  Int_t AnalyzeBuffer(UInt_t* rdata, Bool_t onlysync);
-  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   SetHelicityDetector(THcHelicity* f) { fglHelicityDetector = f; }
-  virtual Int_t  GetNevents() { return fNevents; }
-  virtual Int_t  GetNcycles() { return fNcycles; }
-  virtual Int_t  GetEvNum() { return evNumber; }
-  virtual Int_t* GetHelicityHistoryP() { return fHelicityHistory; }
+  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;}
 
 private:
-  static size_t FindNoCase(const std::string& sdata, const std::string& skey);
 
-  Int_t     fNumBCMs;
-  Double_t* fBCM_Gain;
-  Double_t* fBCM_Offset;
-  Double_t* fBCM_delta_charge;
+  Int_t RanBit30(Int_t ranseed);
+  void MakeParms();
 
-  Int_t  fROC;
   UInt_t fBankID;
   // Helicity Scaler variables
-  Int_t fNevents; /* # of helicity scaler reads in last event */
-  Int_t fNcycles;
+  Int_t fNTrigsInBuf;		/* # of helicity scaler reads in last event */
+  Int_t fNTriggers;
+  Int_t fFirstCycle;
   Int_t fHelicityHistory[200];
   //
   Bool_t fUseFirstEvent;
-  Bool_t fOnlySyncEvents;
-  Bool_t fOnlyBanks;
-  Int_t  fDelayedType;
-  Int_t  fClockChan;
-  UInt_t fLastClock;
-  Int_t  fClockOverflows;
+  Int_t fDelayedType;
+
+  Int_t fRingSeed_reported;
+  Int_t fRingSeed_actual;
+  Int_t fNBits;
+
+  Int_t fNTriggersPlus;
+  Int_t fNTriggersMinus;
+  Double_t* fHScalers[2];
+  Int_t fGateCount[2];
+  Double_t *fScalerSums;
+  Double_t *fAsymmetry;
+  Double_t *fAsymmetryError;
+  Double_t *fCharge;
+  Double_t *fChargeAsymmetry;
+  Double_t fTime;
+  Double_t fTimeAsymmetry;
+  Double_t fTriggerAsymmetry;
 
   std::vector<UInt_t*> fDelayedEvents;
-  std::set<UInt_t>     fRocSet;
+  Int_t fROC;
+  Int_t fNScalerChannels;	// Number of scaler channels/event
+
+  Int_t fNumBCMs;
+  Double_t *fBCM_Gain;
+  Double_t *fBCM_Offset;
+  std::vector <std::string> fBCM_Name;
 
   THcHelicityScaler(const THcHelicityScaler& fh);
   THcHelicityScaler& operator=(const THcHelicityScaler& fh);
 
-  UInt_t       evcount;
-  Double_t     evcountR;
-  UInt_t       evNumber;
-  Int_t        ifound;
-  THcHelicity* fglHelicityDetector;
+  UInt_t evNumber;
+
+  ClassDef(THcHelicityScaler,0)  // Scaler Event handler
 
-  ClassDef(THcHelicityScaler, 0) // Scaler Event handler
 };
 
 #endif
-- 
GitLab