Added dedicated helicity-scaler replays and integrates average BCM currents into standard production replays

CALIBRATION/bcm_current_map/README

@@ -1,19 +1,101 @@
-This script runs over a scaler tree (TSH) and save bcm current values
-for BCM1, BCM2, (BCM4a, BCM4b, BCM17) as well as the event number 
-for each scaler event.
+This script runs over an HMS scaler tree (TSH) or SHMS scaler tree (TSP) and
+saves BCM current values for BCM1, BCM2, BCM4A, BCM4B, and BCM4C, as well as
+the event number for each scaler event.
 
 How to make a parameter file for the average BCM current information
 
-- Use ROOT version 6.08 to use R__LOAD_LIBRARY
-- Do this once to make shared libraries (or after modifying the marco)
+- Do this once to make shared libraries (or after modifying the macro)
 > root -b
 root [0] .L ScalerCalib.C+
 
 To run a script:
-> roob -b
-root [0] .x run.C("/path/Scaler root output file");
+> root -b
+root [0] .L ScalerCalib.C+
+root [1] .x run.C("/path/to/scaler/root/file", "H_or_P");
 
 It prints out 
       # of scaler read, 
       average bcm current values for each bcm
       corresponding event number for each scaler event
+
+
+******************** Comments by rparvez ********************
+Date: Aug 13, 2026
+
+1. With cdaq ROOT 6.36.04, ACLiC places the library under .root_build_dir.
+The bare library name in R__LOAD_LIBRARY is then not found through ROOT's
+library search path. However, doing:
+> root -b
+root [0] .L ScalerCalib.C+
+loads the library in current scope. So, I decided to comment R__LOAD_LIBRARY
+line and adopt following for running the program:
+
+> root -b
+root [0] .L ScalerCalib.C+
+root [1] .x run.C("/path/to/scaler/root/file", "H_or_P");
+
+2. The BCM current workflow uses the dedicated scaler-only replays to make the
+calibration input, then uses the standard production replays to write average
+BCM currents into physics-event tree T.
+
+a. The scaler-only replay scripts are:
+SCRIPTS/COIN/PRODUCTION/replay_helicity_and_scalers_coin.C
+SCRIPTS/HMS/PRODUCTION/replay_helicity_and_scalers_hms.C
+SCRIPTS/SHMS/PRODUCTION/replay_helicity_and_scalers_shms.C
+
+HEEP and HEE runs use these same spectrometer-based scaler replays:
+- HEEP hElec_pProt and pElec_hProt use the COIN replay because both
+  spectrometer arms are read out.
+- HEE hElec uses the HMS replay.
+- HEE pElec uses the SHMS replay.
+
+The particle orientation does not change ordinary or helicity-scaler decoding.
+
+b. Run the appropriate scaler-only replay. Its ROOT filename pattern is:
+COIN: scaler_helicity_replay_coin_<run>_<events>.root
+HMSDIS: scaler_helicity_replay_hms_<run>_<events>.root
+SHMSDIS: scaler_helicity_replay_shms_<run>_<events>.root
+
+c. Run ScalerCalib on the scaler ROOT file to generate bcmcurrent_<run>.param at:
+PARAM/<HMS, SHMS>/BCM/CALIB/bcmcurrent_<run>.param
+note:
+1. I edited ScalerCalib to check and create these subdirectories BCM/, CALIB/.
+2. I decided to use spec_name = "P" for COIN cases, so .param for COIN will go
+to PARAM/SHMS/BCM/CALIB/.
+
+Run ScalerCalib from CALIBRATION/bcm_current_map using:
+
+COIN (reads TSP and writes under PARAM/SHMS):
+root [0] .L ScalerCalib.C+
+root [1] .x run.C("../../ROOTfiles/scaler_helicity_replay_coin_<run>_<events>.root", "P");
+
+HMSDIS (reads TSH and writes under PARAM/HMS):
+root [0] .L ScalerCalib.C+
+root [1] .x run.C("../../ROOTfiles/scaler_helicity_replay_hms_<run>_<events>.root", "H");
+
+SHMSDIS (reads TSP and writes under PARAM/SHMS):
+root [0] .L ScalerCalib.C+
+root [1] .x run.C("../../ROOTfiles/scaler_helicity_replay_shms_<run>_<events>.root", "P");
+
+d. Run the corresponding standard production replay from the repository root:
+
+COIN hElec_pProt:
+hcana -q 'SCRIPTS/COIN/PRODUCTION/replay_production_coin_hElec_pProt.C(<run>,<events>)'
+
+COIN pElec_hProt:
+hcana -q 'SCRIPTS/COIN/PRODUCTION/replay_production_coin_pElec_hProt.C(<run>,<events>)'
+
+HMSDIS or HEE hElec:
+hcana -q 'SCRIPTS/HMS/PRODUCTION/replay_production_hms_coin.C(<run>,<events>)'
+
+SHMSDIS or HEE pElec:
+hcana -q 'SCRIPTS/SHMS/PRODUCTION/replay_production_shms_coin.C(<run>,<events>)'
+
+SetPrintFlag(1) in run.C writes all five BCM current arrays required by
+THcBCMCurrent. The run-specific parameter file is required before starting a
+production replay. The production replay loads that file, maps each scaler
+interval onto physics events, and writes H.bcm.* or P.bcm.* branches in tree T.
+
+If the run-specific parameter file is absent or invalid, the production replay
+stops with a configuration error. Generate and validate the parameter file with
+ScalerCalib before running the corresponding production replay.

CALIBRATION/bcm_current_map/ScalerCalib.C

@@ -1,4 +1,5 @@
 #include <TFile.h>
+#include <TSystem.h>
 #include <TTree.h>
 
 #include "ScalerCalib.h"
@@ -40,21 +41,57 @@ int ScalerCalib::Run()
       return -1;
     }
 
-  int pos;
-  pos = filename.find("scalers_");
-
-  if(pos == -1)
+  string run_token;
+  int pos = filename.find("scaler_helicity_replay_hms_");
+  if(pos != -1)
+    run_token = "scaler_helicity_replay_hms_";
+  else if((pos = filename.find("scaler_helicity_replay_shms_")) != -1)
+    run_token = "scaler_helicity_replay_shms_";
+  else if((pos = filename.find("scaler_helicity_replay_coin_")) != -1)
+    run_token = "scaler_helicity_replay_coin_";
+  else if((pos = filename.find("scaler_helicity_replay_")) != -1)
+    run_token = "scaler_helicity_replay_";
+  else if((pos = filename.find("scalers_")) != -1)
+    run_token = "scalers_";
+  else if((pos = filename.find("production_")) != -1)
+    run_token = "production_";
+  else
     {
-      pos = filename.find("production_");
-      runstr = (filename.substr(pos+11)).substr(0,4);
+      cout << "ERROR: could not determine run number from input filename: "
+           << filename << endl;
+      return -1;
     }
+
+  runstr = filename.substr(pos + run_token.size(), 5);
+
+  string output_dir;
+  if(fName == "H")
+    output_dir = "../../PARAM/HMS/BCM/CALIB";
+  else if(fName == "P")
+    output_dir = "../../PARAM/SHMS/BCM/CALIB";
   else
     {
-      runstr = (filename.substr(pos+8)).substr(0,4);  
+      cout << "ERROR: unsupported spectrometer name: " << fName
+           << ". Use H or P." << endl;
+      return -1;
     }
 
-  ofilename = "bcmcurrent_" + runstr + ".param";
+  if(gSystem->AccessPathName(output_dir.c_str()) &&
+     gSystem->mkdir(output_dir.c_str(), kTRUE) != 0)
+    {
+      cout << "ERROR: could not create output directory: "
+           << output_dir << endl;
+      return -1;
+    }
+
+  ofilename = output_dir + "/bcmcurrent_" + runstr + ".param";
   outfile.open(ofilename.c_str());
+  if(!outfile.is_open())
+    {
+      cout << "ERROR: could not open output parameter file: "
+           << ofilename << endl;
+      return -1;
+    }
 
   outfile << "num_scal_reads = " << evnum.size();
   outfile << "\n" << "\n";

CALIBRATION/bcm_current_map/run.C

@@ -1,4 +1,8 @@
-R__LOAD_LIBRARY(ScalerCalib_C)
+// R__LOAD_LIBRARY(ScalerCalib_C)
+// cdaq ROOT places the ACLiC library under .root_build_dir, while the bare
+// library name (ScalerCalib_C) in run.C is not found through ROOT's search
+// path. So, I commented out this line and load the library from root prompt
+// by running: .L ScalerCalib.C+
 
 void run(string fin="fin.root", string spec_name="H")
 {
@@ -12,4 +16,3 @@ void run(string fin="fin.root", string spec_name="H")
 }
 
 
-

DEF-files/COIN/PRODUCTION/coin_production_hElec_pProt_light.def

@@ -66,6 +66,7 @@ block H.cal.*good*AdcPulseInt
 # *-----------SHMS-------------*
 # ******************************
 # Small blocks
+block P.bcm.*
 block P.ngcer.*
 block P.hgcer.*
 block P.aero.*

DEF-files/COIN/PRODUCTION/coin_production_pElec_hProt_light.def

@@ -66,6 +66,7 @@ block H.cal.*good*AdcPulseInt
 # *-----------SHMS-------------*
 # ******************************
 # Small blocks
+block P.bcm.*
 block P.ngcer.*
 block P.hgcer.*
 block P.aero.*

DEF-files/HMS/EPICS/epics_short.def

@@ -1,4 +1,4 @@
-block T.hms.*
+#block T.hms.*
 
 begin epics
 IBC3H00CRCUR4

DEF-files/HMS/PRODUCTION/BLOCK/hblock_vars_light.def

@@ -16,6 +16,7 @@ block H.kin.*
 block H.rb.*
 block H.react.*
 block H.extcor.*
+block H.bcm.*
 
 # ------------------------------
 # hhodo light

DEF-files/SHMS/EPICS/epics_short.def

@@ -1,4 +1,4 @@
-block T.shms.*
+#block T.shms.*
 
 begin epics
 IBC3H00CRCUR4

DEF-files/SHMS/PRODUCTION/BLOCK/pblock_vars_light.def

@@ -18,6 +18,7 @@ block P.gtr.*
 block P.kin.*
 block P.react.*
 block P.rb.*
+block P.bcm.*
 
 # Commented out prior to light.def
 # block P.tr.*

SCRIPTS/COIN/PRODUCTION/replay_helicity_and_scalers_coin.C

@@ -0,0 +1,146 @@
+#include "MultiFileRun.h"
+
+void replay_helicity_and_scalers_coin(Int_t RunNumber=0, Int_t MaxEvent=0,
+                                      Int_t FirstEvent=1, Int_t MaxSegment=-1) {
+
+  if(RunNumber == 0) {
+    cout << "Enter a Run Number (-1 to exit): ";
+    cin >> RunNumber;
+    if(RunNumber <= 0) return;
+  }
+  if(MaxEvent == 0) {
+    cout << "\nNumber of Events to analyze: ";
+    cin >> MaxEvent;
+    if(MaxEvent == 0) {
+      cerr << "...Invalid entry\n";
+      return;
+    }
+  }
+
+  const char* RunFileNamePattern = "rsidis_production_%05d.dat.%u";
+  vector<string> pathList;
+  pathList.push_back(".");
+  pathList.push_back("./raw");
+  pathList.push_back("./raw/../raw.copiedtotape");
+  pathList.push_back("./cache");
+
+  const char* ROOTFileNamePattern =
+    "ROOTfiles/scaler_helicity_replay_coin_%d_%d.root";
+
+  // Load only the global parameters needed by the scaler handlers.
+  gHcParms->Define("gen_run_number", "Run Number", RunNumber);
+  gHcParms->AddString("g_ctp_database_filename", "DBASE/COIN/standard.database");
+  gHcParms->Load(gHcParms->GetString("g_ctp_database_filename"), RunNumber);
+  gHcParms->Load(gHcParms->GetString("g_ctp_parm_filename"));
+  gHcParms->Load(gHcParms->GetString("g_ctp_kinematics_filename"), RunNumber);
+
+  THcConfigEvtHandler* ev125 =
+    new THcConfigEvtHandler("HC", "Config Event type 125");
+  gHaEvtHandlers->Add(ev125);
+
+  THaEpicsEvtHandler* hcepics =
+    new THaEpicsEvtHandler("epics", "HC EPICS event type 182");
+  gHaEvtHandlers->Add(hcepics);
+
+  THcScalerEvtHandler* pscaler =
+    new THcScalerEvtHandler("P", "Hall C scaler event type 1");
+  pscaler->AddEvtType(1);
+  pscaler->AddEvtType(2);
+  pscaler->AddEvtType(3);
+  pscaler->AddEvtType(4);
+  pscaler->AddEvtType(5);
+  pscaler->AddEvtType(6);
+  pscaler->AddEvtType(7);
+  pscaler->AddEvtType(129);
+  pscaler->SetDelayedType(129);
+  pscaler->SetUseFirstEvent(kTRUE);
+  gHaEvtHandlers->Add(pscaler);
+
+  THcHelicityScaler* phelscaler =
+    new THcHelicityScaler("P", "Hall C helicity scaler");
+  phelscaler->SetROC(8);
+  phelscaler->SetUseFirstEvent(kTRUE);
+  gHaEvtHandlers->Add(phelscaler);
+
+  THcScalerEvtHandler* hscaler =
+    new THcScalerEvtHandler("H", "Hall C scaler event type 2");
+  hscaler->AddEvtType(1);
+  hscaler->AddEvtType(2);
+  hscaler->AddEvtType(3);
+  hscaler->AddEvtType(4);
+  hscaler->AddEvtType(5);
+  hscaler->AddEvtType(6);
+  hscaler->AddEvtType(7);
+  hscaler->AddEvtType(131);
+  hscaler->SetDelayedType(131);
+  hscaler->SetUseFirstEvent(kTRUE);
+  gHaEvtHandlers->Add(hscaler);
+
+  THcHelicityScaler* hhelscaler =
+    new THcHelicityScaler("H", "Hall C helicity scaler");
+  hhelscaler->SetROC(5);
+  hhelscaler->SetUseFirstEvent(kTRUE);
+  gHaEvtHandlers->Add(hhelscaler);
+
+  THcConfigEvtHandler* hconfig =
+    new THcConfigEvtHandler("hconfig", "Hall C configuration event handler");
+  gHaEvtHandlers->Add(hconfig);
+  THcConfigEvtHandler* pconfig =
+    new THcConfigEvtHandler("pconfig", "Hall C configuration event handler");
+  gHaEvtHandlers->Add(pconfig);
+
+  THcAnalyzer* analyzer = new THcAnalyzer;
+  THaEvent* event = new THaEvent;
+
+  vector<string> fileNames = {};
+  Int_t iseg = 0;
+  while(MaxSegment < 0 || iseg <= MaxSegment) {
+    TString codafilename;
+    codafilename.Form(RunFileNamePattern, RunNumber, iseg);
+
+    Bool_t foundSegment = kFALSE;
+    for(UInt_t ipath = 0; ipath < pathList.size(); ipath++) {
+      TString fullPath = TString(pathList[ipath]) + "/" + codafilename;
+      if(!gSystem->AccessPathName(fullPath)) {
+        foundSegment = kTRUE;
+        break;
+      }
+    }
+
+    if(!foundSegment) {
+      if(iseg == 0) {
+        cerr << "ERROR: Required segment 0 file " << codafilename
+             << " was not found in the replay path list. "
+             << "Segment 0 is required to initialize global counters." << endl;
+        cerr << "Searched paths:" << endl;
+        for(UInt_t ipath = 0; ipath < pathList.size(); ipath++)
+          cerr << "  " << pathList[ipath] << endl;
+        return;
+      }
+      cout << "Segment " << iseg << " file " << codafilename
+           << " not found. Finished building input file list." << endl;
+      break;
+    }
+
+    cout << "codafilename = " << codafilename << endl;
+    fileNames.emplace_back(codafilename.Data());
+    iseg++;
+  }
+
+  auto* run = new Podd::MultiFileRun(pathList, fileNames);
+  run->SetRunParamClass("THcRunParameters");
+  run->SetEventRange(FirstEvent, MaxEvent);
+  run->SetNscan(1);
+  run->SetDataRequired(0x7);
+  run->Print();
+
+  TString ROOTFileName = Form(ROOTFileNamePattern, RunNumber, MaxEvent);
+  analyzer->SetCountMode(2);
+  analyzer->SetEvent(event);
+  analyzer->SetEpicsEvtType(182);
+  analyzer->SetCrateMapFileName("MAPS/db_cratemap.dat");
+  analyzer->SetOutFile(ROOTFileName.Data());
+  analyzer->SetOdefFile("DEF-files/COIN/EPICS/epics_short.def");
+  analyzer->EnablePhysicsEvents(false);
+  analyzer->Process(run);
+}