THcHodoscope.cxx

/** \class THcHodoscope
    \ingroup Detectors

\brief Generic hodoscope consisting of multiple
planes with multiple paddles with phototubes on both ends.

This differs from Hall A scintillator class in that it is the whole
hodoscope array, not just one plane.

*/

#include "TClass.h"
#include "THaSubDetector.h"
#include "THcCherenkov.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcRawShowerHit.h"
#include "THcScintPlaneCluster.h"
#include "THcShower.h"
#include "THcSignalHit.h"
#include "math.h"

#include "TClonesArray.h"
#include "THaCutList.h"
#include "THaDetMap.h"
#include "THaEvData.h"
#include "THaGlobals.h"
#include "THaTrack.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THcHodoscope.h"
#include "THcParmList.h"
#include "TMath.h"
#include "VarDef.h"
#include "VarType.h"

#include "THaOutput.h"
#include "TTree.h"

#include "THaTrackProj.h"
#include <vector>

#include <algorithm>
#include <array>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <fstream>
#include <iomanip>
#include <iostream>

#include "hcana/helpers.hxx"
#include <cassert>

using namespace std;
using std::vector;

//_____________________________________________________________________________
THcHodoscope::THcHodoscope(const char* name, const char* description, THaApparatus* apparatus)
    : hcana::ConfigLogging<THaNonTrackingDetector>(name, description, apparatus) {
  // Constructor

  // fTrackProj = new TClonesArray( "THaTrackProj", 5 );
  // Construct the planes
  fNPlanes       = 0; // No planes until we make them
  fStartTime     = -1e5;
  fGoodStartTime = kFALSE;
}

//_____________________________________________________________________________
THcHodoscope::THcHodoscope() : hcana::ConfigLogging<THaNonTrackingDetector>() {
  // Constructor
}

//_____________________________________________________________________________
void THcHodoscope::Setup(const char* name, const char* description) {
  /**
     Create the scintillator plane objects for the hodoscope.

     Uses the Xhodo_num_planes and Xhodo_plane_names to get the number of
     planes and their names.

     Gets a pointer to the Cherenkov named "cer" ("hgcer" in the case of the SHMS.)

  */
  if (IsZombie())
    return;

  // fDebug = 1;  // Keep this at one while we're working on the code

  char prefix[2];

  prefix[0] = tolower(GetApparatus()->GetName()[0]);
  prefix[1] = '\0';

  TString temp(prefix[0]);
  fSHMS = kFALSE;
  if (temp == "p")
    fSHMS = kTRUE;
  TString histname = temp + "_timehist";
  hTime            = new TH1F(histname, "", 400, 0, 200);
  // cout << " fSHMS = " << fSHMS << endl;
  string    planenamelist;
  DBRequest listextra[] = {{"hodo_num_planes", &fNPlanes, kInt},
                           {"hodo_plane_names", &planenamelist, kString},
                           {"hodo_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
                           {"hodo_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
                           {0}};
  // fNPlanes = 4; 		// Default if not defined
  fTDC_RefTimeCut = 0; // Minimum allowed reference times
  fADC_RefTimeCut = 0;
  gHcParms->LoadParmValues((DBRequest*)&listextra, prefix);

  _det_logger->info("Plane Name List : {}", planenamelist);
  // cout << "Plane Name List : " << planenamelist << endl;

  vector<string> plane_names = vsplit(planenamelist);
  // Plane names
  if (plane_names.size() != (UInt_t)fNPlanes) {
    cout << "ERROR: Number of planes " << fNPlanes << " doesn't agree with number of plane names "
         << plane_names.size() << endl;
    // Should quit.  Is there an official way to quit?
  }
  fPlaneNames = new char*[fNPlanes];
  for (Int_t i = 0; i < fNPlanes; i++) {
    fPlaneNames[i] = new char[plane_names[i].length() + 1];
    strcpy(fPlaneNames[i], plane_names[i].c_str());
  }

  // Probably shouldn't assume that description is defined
  char* desc = new char[strlen(description) + 100];
  fPlanes    = new THcScintillatorPlane*[fNPlanes];
  for (Int_t i = 0; i < fNPlanes; i++) {
    strcpy(desc, description);
    strcat(desc, " Plane ");
    strcat(desc, fPlaneNames[i]);
    fPlanes[i] = new THcScintillatorPlane(fPlaneNames[i], desc, i + 1,
                                          this); // Number planes starting from zero!!
    // cout << "Created Scintillator Plane " << fPlaneNames[i] << ", " << desc << endl;
  }

  // Save the nominal particle mass
  THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
  fPartMass                 = app->GetParticleMass();
  fBetaNominal              = app->GetBetaAtPcentral();

  if (fSHMS) {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("hgcer"));
  } else {
    fCherenkov = dynamic_cast<THcCherenkov*>(app->GetDetector("cer"));
  }

  delete[] desc;
}

//_____________________________________________________________________________
THaAnalysisObject::EStatus THcHodoscope::Init(const TDatime& date) {
  // cout << "In THcHodoscope::Init()" << endl;
  Setup(GetName(), GetTitle());

  char EngineDID[] = "xSCIN";
  EngineDID[0]     = toupper(GetApparatus()->GetName()[0]);
  if (gHcDetectorMap->FillMap(fDetMap, EngineDID) < 0) {
    static const char* const here = "Init()";
    // Error( Here(here), "Error filling detectormap for %s.", EngineDID );
    _det_logger->error("THcHodoscope::Init : Error filling detectormap for {}.", EngineDID);
    return kInitError;
  }

  // Should probably put this in ReadDatabase as we will know the
  // maximum number of hits after setting up the detector map
  // But it needs to happen before the sub detectors are initialized
  // so that they can get the pointer to the hitlist.
  _det_logger->info("Hodo TDC and ADC ref time cut = {} {}", fTDC_RefTimeCut, fADC_RefTimeCut);
  // cout << " Hodo tdc ref time cut = " << fTDC_RefTimeCut << " " << fADC_RefTimeCut << endl;

  InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan() + 1, fTDC_RefTimeCut,
              fADC_RefTimeCut);

  EStatus status;
  // This triggers call of ReadDatabase and DefineVariables
  if ((status = THaNonTrackingDetector::Init(date)))
    return fStatus = status;

  for (Int_t ip = 0; ip < fNPlanes; ip++) {
    if ((status = fPlanes[ip]->Init(date))) {
      return fStatus = status;
    }
  }

  // Why not std::vector?
  fNScinHits     = new Int_t[fNPlanes];
  fGoodPlaneTime = new Bool_t[fNPlanes];
  fNPlaneTime    = new Int_t[fNPlanes];
  fSumPlaneTime  = new Double_t[fNPlanes];

  //  Double_t  fHitCnt4 = 0., fHitCnt3 = 0.;

  // Int_t m = 0;
  // fScinHit = new Double_t*[fNPlanes];
  // for ( m = 0; m < fNPlanes; m++ ){
  //   fScinHit[m] = new Double_t[fNPaddle[0]];
  // }

  for (int ip = 0; ip < fNPlanes; ++ip) {
    fScinHitPaddle.push_back(std::vector<Int_t>(fNPaddle[ip], 0));
  }

  fPresentP        = 0;
  THaVar* vpresent = gHaVars->Find(Form("%s.present", GetApparatus()->GetName()));
  if (vpresent) {
    fPresentP = (Bool_t*)vpresent->GetValuePointer();
  }

  return fStatus = kOK;
}
//_____________________________________________________________________________
Int_t THcHodoscope::ReadDatabase(const TDatime& date) {
  /**
     Read this detector's parameters from the ThcParmlist.

     This function is called by THaDetectorBase::Init() once at the
     beginning of the analysis.
  */
  //  static const char* const here = "ReadDatabase()";
  char prefix[2];
  char parname[100];

  // Determine which spectrometer in order to construct
  // the parameter names (e.g. hscin_1x_nr vs. sscin_1x_nr)

  prefix[0] = tolower(GetApparatus()->GetName()[0]);
  //
  prefix[1] = '\0';
  strcpy(parname, prefix);
  strcat(parname, "scin_");
  //  Int_t plen=strlen(parname);
  // cout << " readdatabse hodo fnplanes = " << fNPlanes << endl;

  CreateMissReportParms(Form("%sscin", prefix));

  fBetaNoTrk      = 0.;
  fBetaNoTrkChiSq = 0.;

  fNPaddle      = new UInt_t[fNPlanes];
  fFPTime       = new Double_t[fNPlanes];
  fPlaneCenter  = new Double_t[fNPlanes];
  fPlaneSpacing = new Double_t[fNPlanes];

  prefix[0] = tolower(GetApparatus()->GetName()[0]);
  //
  prefix[1] = '\0';

  for (Int_t i = 0; i < fNPlanes; i++) {

    DBRequest list[] = {{Form("scin_%s_nr", fPlaneNames[i]), &fNPaddle[i], kInt}, {0}};

    gHcParms->LoadParmValues((DBRequest*)&list, prefix);
  }

  // GN added
  // reading variables from *hodo.param
  fMaxScinPerPlane = fNPaddle[0];
  for (Int_t i = 1; i < fNPlanes; i++) {
    fMaxScinPerPlane = (fMaxScinPerPlane > fNPaddle[i]) ? fMaxScinPerPlane : fNPaddle[i];
  }
  // need this for "padded arrays" i.e. 4x16 lists of parameters (GN)
  fMaxHodoScin = fMaxScinPerPlane * fNPlanes;
  if (fDebug >= 1)
    cout << "fMaxScinPerPlane = " << fMaxScinPerPlane << " fMaxHodoScin = " << fMaxHodoScin << endl;

  fHodoVelLight        = new Double_t[fMaxHodoScin];
  fHodoPosSigma        = new Double_t[fMaxHodoScin];
  fHodoNegSigma        = new Double_t[fMaxHodoScin];
  fHodoPosMinPh        = new Double_t[fMaxHodoScin];
  fHodoNegMinPh        = new Double_t[fMaxHodoScin];
  fHodoPosPhcCoeff     = new Double_t[fMaxHodoScin];
  fHodoNegPhcCoeff     = new Double_t[fMaxHodoScin];
  fHodoPosTimeOffset   = new Double_t[fMaxHodoScin];
  fHodoNegTimeOffset   = new Double_t[fMaxHodoScin];
  fHodoPosPedLimit     = new Int_t[fMaxHodoScin];
  fHodoNegPedLimit     = new Int_t[fMaxHodoScin];
  fHodoPosInvAdcOffset = new Double_t[fMaxHodoScin];
  fHodoNegInvAdcOffset = new Double_t[fMaxHodoScin];
  fHodoPosInvAdcLinear = new Double_t[fMaxHodoScin];
  fHodoNegInvAdcLinear = new Double_t[fMaxHodoScin];
  fHodoPosInvAdcAdc    = new Double_t[fMaxHodoScin];
  fHodoNegInvAdcAdc    = new Double_t[fMaxHodoScin];

  // New Time-Walk Calibration Parameters
  fHodoVelFit   = new Double_t[fMaxHodoScin];
  fHodoCableFit = new Double_t[fMaxHodoScin];
  fHodo_LCoeff  = new Double_t[fMaxHodoScin];
  fHodoPos_c1   = new Double_t[fMaxHodoScin];
  fHodoNeg_c1   = new Double_t[fMaxHodoScin];
  fHodoPos_c2   = new Double_t[fMaxHodoScin];
  fHodoNeg_c2   = new Double_t[fMaxHodoScin];
  fHodoSigmaPos = new Double_t[fMaxHodoScin];
  fHodoSigmaNeg = new Double_t[fMaxHodoScin];

  fNHodoscopes             = 2;
  fxLoScin                 = new Int_t[fNHodoscopes];
  fxHiScin                 = new Int_t[fNHodoscopes];
  fyLoScin                 = new Int_t[fNHodoscopes];
  fyHiScin                 = new Int_t[fNHodoscopes];
  fHodoSlop                = new Double_t[fNPlanes];
  fTdcOffset               = new Int_t[fNPlanes];
  fAdcTdcOffset            = new Double_t[fNPlanes];
  fHodoPosAdcTimeWindowMin = new Double_t[fMaxHodoScin];
  fHodoPosAdcTimeWindowMax = new Double_t[fMaxHodoScin];
  fHodoNegAdcTimeWindowMin = new Double_t[fMaxHodoScin];
  fHodoNegAdcTimeWindowMax = new Double_t[fMaxHodoScin];

  for (Int_t ip = 0; ip < fNPlanes; ip++) { // Set a large default window
    fTdcOffset[ip]    = 0;
    fAdcTdcOffset[ip] = 0.0;
  }

  DBRequest list[] = {
      {"cosmicflag", &fCosmicFlag, kInt, 0, 1},
      {"NumPlanesBetaCalc", &fNumPlanesBetaCalc, kInt, 0, 1},
      {"start_time_center", &fStartTimeCenter, kDouble},
      {"start_time_slop", &fStartTimeSlop, kDouble},
      {"scin_tdc_to_time", &fScinTdcToTime, kDouble},
      {"scin_tdc_min", &fScinTdcMin, kDouble},
      {"scin_tdc_max", &fScinTdcMax, kDouble},
      {"tof_tolerance", &fTofTolerance, kDouble, 0, 1},
      {"pathlength_central", &fPathLengthCentral, kDouble},
      {"hodo_pos_sigma", &fHodoPosSigma[0], kDouble, fMaxHodoScin, 1},
      {"hodo_neg_sigma", &fHodoNegSigma[0], kDouble, fMaxHodoScin, 1},
      {"hodo_pos_ped_limit", &fHodoPosPedLimit[0], kInt, fMaxHodoScin, 1},
      {"hodo_neg_ped_limit", &fHodoNegPedLimit[0], kInt, fMaxHodoScin, 1},
      {"tofusinginvadc", &fTofUsingInvAdc, kInt, 0, 1},
      {"xloscin", &fxLoScin[0], kInt, (UInt_t)fNHodoscopes},
      {"xhiscin", &fxHiScin[0], kInt, (UInt_t)fNHodoscopes},
      {"yloscin", &fyLoScin[0], kInt, (UInt_t)fNHodoscopes},
      {"yhiscin", &fyHiScin[0], kInt, (UInt_t)fNHodoscopes},
      {"track_eff_test_num_scin_planes", &fTrackEffTestNScinPlanes, kInt},
      {"trackeff_scint_ydiff_max", &trackeff_scint_ydiff_max, kDouble, 0, 1},
      {"trackeff_scint_xdiff_max", &trackeff_scint_xdiff_max, kDouble, 0, 1},
      {"cer_npe", &fNCerNPE, kDouble, 0, 1},
      {"normalized_energy_tot", &fNormETot, kDouble, 0, 1},
      {"hodo_slop", fHodoSlop, kDouble, (UInt_t)fNPlanes},
      {"debugprintscinraw", &fdebugprintscinraw, kInt, 0, 1},
      {"hodo_tdc_offset", fTdcOffset, kInt, (UInt_t)fNPlanes, 1},
      {"hodo_adc_tdc_offset", fAdcTdcOffset, kDouble, (UInt_t)fNPlanes, 1},
      {"hodo_PosAdcTimeWindowMin", fHodoPosAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
      {"hodo_PosAdcTimeWindowMax", fHodoPosAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
      {"hodo_NegAdcTimeWindowMin", fHodoNegAdcTimeWindowMin, kDouble, (UInt_t)fMaxHodoScin, 1},
      {"hodo_NegAdcTimeWindowMax", fHodoNegAdcTimeWindowMax, kDouble, (UInt_t)fMaxHodoScin, 1},
      {"dumptof", &fDumpTOF, kInt, 0, 1},
      {"TOFCalib_shtrk_lo", &fTOFCalib_shtrk_lo, kDouble, 0, 1},
      {"TOFCalib_shtrk_hi", &fTOFCalib_shtrk_hi, kDouble, 0, 1},
      {"TOFCalib_cer_lo", &fTOFCalib_cer_lo, kDouble, 0, 1},
      {"TOFCalib_beta_lo", &fTOFCalib_beta_lo, kDouble, 0, 1},
      {"TOFCalib_beta_hi", &fTOFCalib_beta_hi, kDouble, 0, 1},
      {"dumptof_filename", &fTOFDumpFile, kString, 0, 1},
      {0}};

  // Defaults if not defined in parameter file

  trackeff_scint_ydiff_max = 20.;
  trackeff_scint_xdiff_max = 20.;
  for (UInt_t ip = 0; ip < fMaxHodoScin; ip++) {
    fHodoPosAdcTimeWindowMin[ip] = -1000.;
    fHodoPosAdcTimeWindowMax[ip] = 1000.;
    fHodoNegAdcTimeWindowMin[ip] = -1000.;
    fHodoNegAdcTimeWindowMax[ip] = 1000.;

    fHodoPosPedLimit[ip] = 0.0;
    fHodoNegPedLimit[ip] = 0.0;
    fHodoPosSigma[ip]    = 0.2;
    fHodoNegSigma[ip]    = 0.2;
  }
  fTOFCalib_shtrk_lo = -kBig;
  fTOFCalib_shtrk_hi = kBig;
  fTOFCalib_cer_lo   = -kBig;
  fTOFCalib_beta_lo  = -kBig;
  fTOFCalib_beta_hi  = kBig;
  fdebugprintscinraw = 0;
  fDumpTOF           = 0;
  fTOFDumpFile       = "";
  fTofUsingInvAdc    = 1;
  fTofTolerance      = 3.0;
  fNCerNPE           = 2.0;
  fNormETot          = 0.7;
  fCosmicFlag        = 0;
  fNumPlanesBetaCalc = 4;
  // Gets added to each reference time corrected raw TDC value
  // to make sure valid range is all positive.

  gHcParms->LoadParmValues((DBRequest*)&list, prefix);
  if (fCosmicFlag == 1)
    cout << "Setup for cosmics in TOF" << endl;
  // cout << " cosmic flag = " << fCosmicFlag << endl;
  if (fDumpTOF) {
    fDumpOut.open(fTOFDumpFile.c_str());
    if (fDumpOut.is_open()) {
      // fDumpOut << "Hodoscope Time of Flight calibration data" << endl;
    } else {
      fDumpTOF = 0;
      cout << "WARNING: Unable to open TOF Dump file " << fTOFDumpFile << endl;
      cout << "Data for TOF calibration not being written." << endl;
    }
  }

  // cout << " x1 lo = " << fxLoScin[0]
  //      << " x2 lo = " << fxLoScin[1]
  //      << " x1 hi = " << fxHiScin[0]
  //      << " x2 hi = " << fxHiScin[1]
  //      << endl;

  // cout << " y1 lo = " << fyLoScin[0]
  //      << " y2 lo = " << fyLoScin[1]
  //      << " y1 hi = " << fyHiScin[0]
  //      << " y2 hi = " << fyHiScin[1]
  //      << endl;

  // cout << "Hdososcope planes hits for trigger = " << fTrackEffTestNScinPlanes
  //      << " normalized energy min = " << fNormETot
  //      << " number of photo electrons = " << fNCerNPE
  //      << endl;

  // Set scin Velocity/Cable to default
  for (UInt_t i = 0; i < fMaxHodoScin; i++) {
    fHodoVelLight[i] = 15.0;
  }

  if (fTofUsingInvAdc) {
    DBRequest list2[] = {
        {"hodo_vel_light", &fHodoVelLight[0], kDouble, fMaxHodoScin, 1},
        {"hodo_pos_invadc_offset", &fHodoPosInvAdcOffset[0], kDouble, fMaxHodoScin},
        {"hodo_neg_invadc_offset", &fHodoNegInvAdcOffset[0], kDouble, fMaxHodoScin},
        {"hodo_pos_invadc_linear", &fHodoPosInvAdcLinear[0], kDouble, fMaxHodoScin},
        {"hodo_neg_invadc_linear", &fHodoNegInvAdcLinear[0], kDouble, fMaxHodoScin},
        {"hodo_pos_invadc_adc", &fHodoPosInvAdcAdc[0], kDouble, fMaxHodoScin},
        {"hodo_neg_invadc_adc", &fHodoNegInvAdcAdc[0], kDouble, fMaxHodoScin},
        {0}};

    gHcParms->LoadParmValues((DBRequest*)&list2, prefix);
  };
  /* if (!fTofUsingInvAdc) {
      DBRequest list3[]={
    {"hodo_vel_light",                   &fHodoVelLight[0],       kDouble,  fMaxHodoScin},
    {"hodo_pos_minph",                   &fHodoPosMinPh[0],       kDouble,  fMaxHodoScin},
    {"hodo_neg_minph",                   &fHodoNegMinPh[0],       kDouble,  fMaxHodoScin},
    {"hodo_pos_phc_coeff",               &fHodoPosPhcCoeff[0],    kDouble,  fMaxHodoScin},
    {"hodo_neg_phc_coeff",               &fHodoNegPhcCoeff[0],    kDouble,  fMaxHodoScin},
    {"hodo_pos_time_offset",             &fHodoPosTimeOffset[0],  kDouble,  fMaxHodoScin},
    {"hodo_neg_time_offset",             &fHodoNegTimeOffset[0],  kDouble,  fMaxHodoScin},
    {0}
    };


    gHcParms->LoadParmValues((DBRequest*)&list3,prefix);
  */
  DBRequest list4[] = {{"hodo_velFit", &fHodoVelFit[0], kDouble, fMaxHodoScin, 1},
                       {"hodo_cableFit", &fHodoCableFit[0], kDouble, fMaxHodoScin, 1},
                       {"hodo_LCoeff", &fHodo_LCoeff[0], kDouble, fMaxHodoScin, 1},
                       {"c1_Pos", &fHodoPos_c1[0], kDouble, fMaxHodoScin, 1},
                       {"c1_Neg", &fHodoNeg_c1[0], kDouble, fMaxHodoScin, 1},
                       {"c2_Pos", &fHodoPos_c2[0], kDouble, fMaxHodoScin, 1},
                       {"c2_Neg", &fHodoNeg_c2[0], kDouble, fMaxHodoScin, 1},
                       {"TDC_threshold", &fTdc_Thrs, kDouble, 0, 1},
                       {"hodo_PosSigma", &fHodoSigmaPos[0], kDouble, fMaxHodoScin, 1},
                       {"hodo_NegSigma", &fHodoSigmaNeg[0], kDouble, fMaxHodoScin, 1},
                       {0}};

  fTdc_Thrs = 1.0;
  // Set Default Values if NOT defined in param file
  for (UInt_t i = 0; i < fMaxHodoScin; i++) {

    // Turn OFF Time-Walk Correction if param file NOT found
    fHodoPos_c1[i] = 0.0;
    fHodoPos_c2[i] = 0.0;
    fHodoNeg_c1[i] = 0.0;
    fHodoNeg_c2[i] = 0.0;
  }
  for (UInt_t i = 0; i < fMaxHodoScin; i++) {
    // Set scin Velocity/Cable to default
    fHodoCableFit[i] = 0.0;
    fHodoVelFit[i]   = 15.0;
    // set time coeff between paddles to default
    fHodo_LCoeff[i] = 0.0;
  }

  gHcParms->LoadParmValues((DBRequest*)&list4, prefix);

  if (fDebug >= 1) {
    cout << "******* Testing Hodoscope Parameter Reading ***\n";
    cout << "StarTimeCenter = " << fStartTimeCenter << endl;
    cout << "StartTimeSlop = " << fStartTimeSlop << endl;
    cout << "ScintTdcToTime = " << fScinTdcToTime << endl;
    cout << "TdcMin = " << fScinTdcMin << " TdcMax = " << fScinTdcMax << endl;
    cout << "TofTolerance = " << fTofTolerance << endl;
    cout << "*** VelLight ***\n";
    for (Int_t i1 = 0; i1 < fNPlanes; i1++) {
      cout << "Plane " << i1 << endl;
      for (UInt_t i2 = 0; i2 < fMaxScinPerPlane; i2++) {
        cout << fHodoVelLight[GetScinIndex(i1, i2)] << " ";
      }
      cout << endl;
    }
    cout << endl << endl;
    // check fHodoPosPhcCoeff
    /*
    cout <<"fHodoPosPhcCoeff = ";
    for (int i1=0;i1<fMaxHodoScin;i1++) {
      cout<<this->GetHodoPosPhcCoeff(i1)<<" ";
    }
    cout<<endl;
    */
  }
  //
  if ((fTofTolerance > 0.5) && (fTofTolerance < 10000.)) {
    // cout << "USING "<<fTofTolerance<<" NSEC WINDOW FOR FP NO_TRACK CALCULATIONS.\n";
    _det_logger->info("THcHodoscope: Using {} nsec window for fp no_track calculations.",
                      fTofTolerance);
  } else {
    fTofTolerance = 3.0;
    // cout << "*** USING DEFAULT 3 NSEC WINDOW FOR FP NO_TRACK CALCULATIONS!! ***\n";
    _det_logger->warn("THcHodoscope: Using default {} nsec window for fp no_track calculations.",
                      fTofTolerance);
  }
  fRatio_xpfp_to_xfp = 0.00;
  TString SHMS       = "p";
  TString HMS        = "h";
  TString test       = prefix[0];
  if (test == SHMS)
    fRatio_xpfp_to_xfp = 0.0018; // SHMS
  if (test == HMS)
    fRatio_xpfp_to_xfp = 0.0011; // HMS
  cout << " fRatio_xpfp_to_xfp= " << fRatio_xpfp_to_xfp << endl;

  fIsInit = true;
  return kOK;
}

//_____________________________________________________________________________
Int_t THcHodoscope::DefineVariables(EMode mode) {
  /**
    Initialize global variables for histograms and Root tree
  */
  // cout << "THcHodoscope::DefineVariables called " << GetName() << endl;
  if (mode == kDefine && fIsSetup)
    return kOK;
  fIsSetup = (mode == kDefine);

  // Register variables in global list

  RVarDef vars[] = {
      // Move these into THcHallCSpectrometer using track fTracks
      {"beta", "Beta including track info", "fBeta"},
      {"betanotrack", "Beta from scintillator hits", "fBetaNoTrk"},
      {"betachisqnotrack", "Chi square of beta from scintillator hits", "fBetaNoTrkChiSq"},
      {"fpHitsTime", "Time at focal plane from all hits", "fFPTimeAll"},
      {"starttime", "Hodoscope Start Time", "fStartTime"},
      {"goodstarttime", "Hodoscope Good Start Time (logical flag)", "fGoodStartTime"},
      {"goodscinhit", "Hit in fid area", "fGoodScinHits"},
      {"TimeHist_Sigma", "", "fTimeHist_Sigma"},
      {"TimeHist_Peak", "", "fTimeHist_Peak"},
      {"TimeHist_Hits", "", "fTimeHist_Hits"},
      {0}};

  return DefineVarsFromList(vars, mode);
  //  return kOK;
}
//_____________________________________________________________________________

Int_t THcHodoscope::ManualInitTree(TTree* t) {
  // The most direct path to the output tree!!!
  std::string app_name    = GetApparatus()->GetName();
  std::string det_name    = GetName();
  std::string branch_name = (app_name + "_" + det_name + "_data");
  if (t) {
    _det_logger->info("THcHodoscope::ManualInitTree : Adding branch, {}, to output tree",
                      branch_name);
    t->Branch(branch_name.c_str(), &_basic_data, 32000, 99);
  }
  return 0;
}
//_____________________________________________________________________________

THcHodoscope::~THcHodoscope() {
  // Destructor. Remove variables from global list.

  delete[] fFPTime;
  delete[] fPlaneCenter;
  delete[] fPlaneSpacing;

  if (fIsSetup)
    RemoveVariables();
  if (fIsInit)
    DeleteArrays();

  for (int i = 0; i < fNPlanes; ++i) {
    delete fPlanes[i];
    delete[] fPlaneNames[i];
  }
  delete[] fPlanes;
  delete[] fPlaneNames;
}

//_____________________________________________________________________________
void THcHodoscope::DeleteArrays() {
  // Delete member arrays. Used by destructor.
  // Int_t k;
  // for( k = 0; k < fNPlanes; k++){
  //   delete [] fScinHit[k];
  // }
  // delete [] fScinHit;

  delete[] fxLoScin;
  fxLoScin = NULL;
  delete[] fxHiScin;
  fxHiScin = NULL;
  delete[] fyLoScin;
  fyLoScin = NULL;
  delete[] fyHiScin;
  fyHiScin = NULL;
  delete[] fHodoSlop;
  fHodoSlop = NULL;

  delete[] fNPaddle;
  fNPaddle = NULL;
  delete[] fHodoVelLight;
  fHodoVelLight = NULL;
  delete[] fHodoPosSigma;
  fHodoPosSigma = NULL;
  delete[] fHodoNegSigma;
  fHodoNegSigma = NULL;
  delete[] fHodoPosMinPh;
  fHodoPosMinPh = NULL;
  delete[] fHodoNegMinPh;
  fHodoNegMinPh = NULL;
  delete[] fHodoPosPhcCoeff;
  fHodoPosPhcCoeff = NULL;
  delete[] fHodoNegPhcCoeff;
  fHodoNegPhcCoeff = NULL;
  delete[] fHodoPosTimeOffset;
  fHodoPosTimeOffset = NULL;
  delete[] fHodoNegTimeOffset;
  fHodoNegTimeOffset = NULL;
  delete[] fHodoPosPedLimit;
  fHodoPosPedLimit = NULL;
  delete[] fHodoNegPedLimit;
  fHodoNegPedLimit = NULL;
  delete[] fHodoPosInvAdcOffset;
  fHodoPosInvAdcOffset = NULL;
  delete[] fHodoNegInvAdcOffset;
  fHodoNegInvAdcOffset = NULL;
  delete[] fHodoPosInvAdcLinear;
  fHodoPosInvAdcLinear = NULL;
  delete[] fHodoNegInvAdcLinear;
  fHodoNegInvAdcLinear = NULL;
  delete[] fHodoPosInvAdcAdc;
  fHodoPosInvAdcAdc = NULL;
  delete[] fHodoNegInvAdcAdc;
  fHodoNegInvAdcAdc = NULL;
  delete[] fGoodPlaneTime;
  fGoodPlaneTime = NULL;
  delete[] fNPlaneTime;
  fNPlaneTime = NULL;
  delete[] fSumPlaneTime;
  fSumPlaneTime = NULL;
  delete[] fNScinHits;
  fNScinHits = NULL;
  delete[] fTdcOffset;
  fTdcOffset = NULL;
  delete[] fAdcTdcOffset;
  fAdcTdcOffset = NULL;
  delete[] fHodoNegAdcTimeWindowMin;
  fHodoNegAdcTimeWindowMin = NULL;
  delete[] fHodoNegAdcTimeWindowMax;
  fHodoNegAdcTimeWindowMax = NULL;
  delete[] fHodoPosAdcTimeWindowMin;
  fHodoPosAdcTimeWindowMin = NULL;
  delete[] fHodoPosAdcTimeWindowMax;
  fHodoPosAdcTimeWindowMax = NULL;

  delete[] fHodoVelFit;
  fHodoVelFit = NULL;
  delete[] fHodoCableFit;
  fHodoCableFit = NULL;
  delete[] fHodo_LCoeff;
  fHodo_LCoeff = NULL;
  delete[] fHodoPos_c1;
  fHodoPos_c1 = NULL;
  delete[] fHodoNeg_c1;
  fHodoNeg_c1 = NULL;
  delete[] fHodoPos_c2;
  fHodoPos_c2 = NULL;
  delete[] fHodoNeg_c2;
  fHodoNeg_c2 = NULL;
  delete[] fHodoSigmaPos;
  fHodoSigmaPos = NULL;
  delete[] fHodoSigmaNeg;
  fHodoSigmaNeg = NULL;
}

//_____________________________________________________________________________
void THcHodoscope::Clear(Option_t* opt) {
  /*! \brief Clears variables
   *
   *  Called by  THcHodoscope::Decode
   *
   */
  fTimeHist_Sigma = kBig;
  fTimeHist_Peak  = kBig;
  fTimeHist_Hits  = kBig;

  fBeta           = 0.0;
  fBetaNoTrk      = 0.0;
  fBetaNoTrkChiSq = 0.0;
  fStartTime      = -1000.;
  fFPTimeAll      = -1000.;
  fGoodStartTime  = kFALSE;
  fGoodScinHits   = 0;

  if (*opt != 'I') {
    for (Int_t ip = 0; ip < fNPlanes; ip++) {
      fPlanes[ip]->Clear();
      fFPTime[ip]       = 0.;
      fPlaneCenter[ip]  = 0.;
      fPlaneSpacing[ip] = 0.;
      for (UInt_t iPaddle = 0; iPaddle < fNPaddle[ip]; ++iPaddle) {
        fScinHitPaddle[ip][iPaddle] = 0;
      }
    }
  }
  fdEdX.clear();
  fNScinHit.clear();
  fNClust.clear();
  fClustSize.clear();
  fClustPos.clear();
  fThreeScin.clear();
  fGoodScinHitsX.clear();
  fGoodFlags.clear();
}

//_____________________________________________________________________________
Int_t THcHodoscope::Decode(const THaEvData& evdata) {
  /*! \brief Decodes raw data and processes raw data into hits for each instance of
   * THcScintillatorPlane
   *
   *  - Reads raw data using THcHitList::DecodeToHitList
   *  - If one wants to subtract pedestals (assumed to be a set of data at beginning of run)
   *    + Must define "Pedestal_event" cut in the cuts definition file
   *    + For each "Pedestal_event" calls THcScintillatorPlane::AccumulatePedestals and returns
   *    + After First event which is not a  "Pedestal_event" calls
   * THcScintillatorPlane::CalculatePedestals
   *  - For each scintillator plane THcScintillatorPlane::ProcessHits
   *  - Calls THcHodoscope::EstimateFocalPlaneTime
   *
   *
   */
  // Get the Hall C style hitlist (fRawHitList) for this event
  Bool_t present = kTRUE; // Suppress reference time warnings
  if (fPresentP) {        // if this spectrometer not part of trigger
    present = *fPresentP;
  }
  fNHits = DecodeToHitList(evdata, !present);

  //
  // GN: print event number so we can cross-check with engine
  // if (evdata.GetEvNum()>1000)
  //   cout <<"\nhcana_event " << evdata.GetEvNum()<<endl;

  fCheckEvent = evdata.GetEvNum();
  fEventType  = evdata.GetEvType();

  if (gHaCuts->Result("Pedestal_event")) {
    Int_t nexthit = 0;
    for (Int_t ip = 0; ip < fNPlanes; ip++) {

      nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
    }
    fAnalyzePedestals = 1; // Analyze pedestals first normal events
    return (0);
  }
  if (fAnalyzePedestals) {
    for (Int_t ip = 0; ip < fNPlanes; ip++) {

      fPlanes[ip]->CalculatePedestals();
    }
    fAnalyzePedestals = 0; // Don't analyze pedestals next event
  }

  // Let each plane get its hits
  Int_t nexthit = 0;

  fNfptimes   = 0;
  Int_t thits = 0;
  for (Int_t ip = 0; ip < fNPlanes; ip++) {

    fPlaneCenter[ip]  = fPlanes[ip]->GetPosCenter(0) + fPlanes[ip]->GetPosOffset();
    fPlaneSpacing[ip] = fPlanes[ip]->GetSpacing();

    //    nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
    // GN: select only events that have reasonable TDC values to start with
    // as per the Engine h_strip_scin.f
    nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
    thits += fPlanes[ip]->GetNScinHits();
  }
  fStartTime = -1000;
  if (thits > 0)
    EstimateFocalPlaneTime();

  if (fdebugprintscinraw == 1) {
    for (UInt_t ihit = 0; ihit < fNRawHits; ihit++) {
      //    THcRawHodoHit* hit = (THcRawHodoHit *) fRawHitList->At(ihit);
      //    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
      //	 << hit->fADC_pos << " " << hit->fADC_neg << " "  <<  hit->fTDC_pos
      //	 << " " <<  hit->fTDC_neg << endl;
    }
    cout << endl;
  }

  return fNHits;
}

//_____________________________________________________________________________
void THcHodoscope::EstimateFocalPlaneTime() {
  /*! \brief Calculates the Drift Chamber start time and fBetaNoTrk (velocity determined without
   * track info)
   *
   *  - Called by  THcHodoscope::Decode
   *  - selects good scintillator paddle hits
   *     + loops through hits in each scintillator plane and fills histogram array, "timehist", with
   * corrected times for positive and negative ends of each paddle
   *     + Determines the peak of "timehist"
   *
   */
  Int_t ihit      = 0;
  Int_t nscinhits = 0; // Total # hits with at least one good tdc
  hTime->Reset();
  //
  for (Int_t ip = 0; ip < fNPlanes; ip++) {
    Int_t nphits = fPlanes[ip]->GetNScinHits();
    nscinhits += nphits;
    TClonesArray* hodoHits = fPlanes[ip]->GetHits();
    for (Int_t i = 0; i < nphits; i++) {
      THcHodoHit* hit = (THcHodoHit*)hodoHits->At(i);
      if (hit->GetHasCorrectedTimes()) {
        Double_t postime = hit->GetPosTOFCorrectedTime();
        Double_t negtime = hit->GetNegTOFCorrectedTime();
        hTime->Fill(postime);
        hTime->Fill(negtime);
      }
    }
  }
  //
  //
  ihit                      = 0;
  Double_t fpTimeSum        = 0.0;
  fNfptimes                 = 0;
  Int_t    Ngood_hits_plane = 0;
  Double_t Plane_fptime_sum = 0.0;
  Bool_t   goodplanetime[fNPlanes];
  Bool_t   twogoodtimes[nscinhits];
  Double_t tmin               = 0.5 * hTime->GetMaximumBin();
  fTimeHist_Peak              = tmin;
  fTimeHist_Sigma             = hTime->GetRMS();
  fTimeHist_Hits              = hTime->Integral();
  _basic_data.fTimeHist_Peak  = fTimeHist_Peak;
  _basic_data.fTimeHist_Sigma = fTimeHist_Sigma;
  _basic_data.fTimeHist_Hits  = fTimeHist_Hits;
  for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
    goodplanetime[ip]      = kFALSE;
    Int_t         nphits   = fPlanes[ip]->GetNScinHits();
    TClonesArray* hodoHits = fPlanes[ip]->GetHits();
    Ngood_hits_plane       = 0;
    Plane_fptime_sum       = 0.0;
    for (Int_t i = 0; i < nphits; i++) {
      THcHodoHit* hit    = (THcHodoHit*)hodoHits->At(i);
      twogoodtimes[ihit] = kFALSE;
      if (hit->GetHasCorrectedTimes()) {
        Double_t postime = hit->GetPosTOFCorrectedTime();
        Double_t negtime = hit->GetNegTOFCorrectedTime();
        if ((postime > (tmin - fTofTolerance)) && (postime < (tmin + fTofTolerance)) &&
            (negtime > (tmin - fTofTolerance)) && (negtime < (tmin + fTofTolerance))) {
          hit->SetTwoGoodTimes(kTRUE);
          twogoodtimes[ihit] = kTRUE;                      // Both tubes fired
          Int_t    index     = hit->GetPaddleNumber() - 1; //
          Double_t fptime;
          if (fCosmicFlag == 1) {
            fptime = hit->GetScinCorrectedTime() +
                     (fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
                         (29.979 * fBetaNominal);
          } else {
            fptime = hit->GetScinCorrectedTime() -
                     (fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos()) /
                         (29.979 * fBetaNominal);
          }
          Ngood_hits_plane++;
          Plane_fptime_sum += fptime;
          fpTimeSum += fptime;
          fNfptimes++;
          goodplanetime[ip] = kTRUE;
        } else {
          hit->SetTwoGoodTimes(kFALSE);
        }
      }
      ihit++;
    }
    if (Ngood_hits_plane)
      fPlanes[ip]->SetFpTime(Plane_fptime_sum / float(Ngood_hits_plane));
    fPlanes[ip]->SetNGoodHits(Ngood_hits_plane);
  }

  if (fNfptimes > 0) {
    fStartTime     = fpTimeSum / fNfptimes;
    fGoodStartTime = kTRUE;
    fFPTimeAll     = fStartTime;
  } else {
    fStartTime     = fStartTimeCenter;
    fGoodStartTime = kFALSE;
    fFPTimeAll     = fStartTime;
  }
  //
  //
  hTime->Reset();
  //
  if ((goodplanetime[0] || goodplanetime[1]) && (goodplanetime[2] || goodplanetime[3])) {

    Double_t sumW  = 0.;
    Double_t sumT  = 0.;
    Double_t sumZ  = 0.;
    Double_t sumZZ = 0.;
    Double_t sumTZ = 0.;
    Int_t    ihhit = 0;

    for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) {
      Int_t         nphits   = fPlanes[ip]->GetNScinHits();
      TClonesArray* hodoHits = fPlanes[ip]->GetHits();

      for (Int_t i = 0; i < nphits; i++) {
        Int_t index = ((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber() - 1;

        if (twogoodtimes[ihhit]) {

          Double_t sigma = 0.0;
          if (fTofUsingInvAdc) {
            sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoPosSigma[GetScinIndex(ip, index)], 2) +
                                       TMath::Power(fHodoNegSigma[GetScinIndex(ip, index)], 2)));
          } else {
            sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoSigmaPos[GetScinIndex(ip, index)], 2) +
                                       TMath::Power(fHodoSigmaNeg[GetScinIndex(ip, index)], 2)));
          }

          Double_t scinWeight = 1 / TMath::Power(sigma, 2);
          Double_t zPosition  = fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos();
          //	  cout << "hit = " << ihhit + 1 << "   zpos = " << zPosition << "   sigma = " <<
          // sigma << endl; cout << "fHodoSigma+ = " << fHodoSigmaPos[GetScinIndex(ip,index)] <<
          // endl;
          sumW += scinWeight;
          sumT += scinWeight * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();
          sumZ += scinWeight * zPosition;
          sumZZ += scinWeight * (zPosition * zPosition);
          sumTZ += scinWeight * zPosition * ((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime();

        } // condition of good scin time
        ihhit++;
      } // loop over hits of plane
    }   // loop over planes

    Double_t tmp      = sumW * sumZZ - sumZ * sumZ;
    Double_t t0       = (sumT * sumZZ - sumZ * sumTZ) / tmp;
    Double_t tmpDenom = sumW * sumTZ - sumZ * sumT;

    if (TMath::Abs(tmpDenom) > (1 / 10000000000.0)) {

      fBetaNoTrk      = tmp / tmpDenom;
      fBetaNoTrkChiSq = 0.;
      ihhit           = 0;

      for (Int_t ip = 0; ip < fNumPlanesBetaCalc; ip++) { // Loop over planes
        Int_t         nphits   = fPlanes[ip]->GetNScinHits();
        TClonesArray* hodoHits = fPlanes[ip]->GetHits();

        for (Int_t i = 0; i < nphits; i++) {
          Int_t index = ((THcHodoHit*)hodoHits->At(i))->GetPaddleNumber() - 1;

          if (twogoodtimes[ihhit]) {

            Double_t zPosition = fPlanes[ip]->GetZpos() + (index % 2) * fPlanes[ip]->GetDzpos();
            Double_t timeDif   = (((THcHodoHit*)hodoHits->At(i))->GetScinCorrectedTime() - t0);

            Double_t sigma = 0.0;
            if (fTofUsingInvAdc) {
              sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoPosSigma[GetScinIndex(ip, index)], 2) +
                                         TMath::Power(fHodoNegSigma[GetScinIndex(ip, index)], 2)));
            } else {
              sigma = 0.5 * (TMath::Sqrt(TMath::Power(fHodoSigmaPos[GetScinIndex(ip, index)], 2) +
                                         TMath::Power(fHodoSigmaNeg[GetScinIndex(ip, index)], 2)));
            }

            fBetaNoTrkChiSq +=
                ((zPosition / fBetaNoTrk - timeDif) * (zPosition / fBetaNoTrk - timeDif)) /
                (sigma * sigma);