THcCherenkov.cxx

/** \class THcCherenkov
    \ingroup Detectors

\brief Class for gas Cherenkov detectors

*/

#include "THcCherenkov.h"
#include "TClonesArray.h"
#include "THaApparatus.h"
#include "THaCutList.h"
#include "THaDetMap.h"
#include "THaEvData.h"
#include "THaTrack.h"
#include "THaTrackProj.h"
#include "THcDetectorMap.h"
#include "THcGlobals.h"
#include "THcHallCSpectrometer.h"
#include "THcHitList.h"
#include "THcHodoscope.h"
#include "THcParmList.h"
#include "THcSignalHit.h"
#include "TMath.h"
#include "VarDef.h"
#include "VarType.h"

#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <iomanip>
#include <iostream>
#include <string>

using namespace std;

using std::cin;
using std::cout;
using std::endl;
using std::setprecision;
using std::setw;

//_____________________________________________________________________________
THcCherenkov::THcCherenkov(const char* name, const char* description, THaApparatus* apparatus)
    : THaNonTrackingDetector(name, description, apparatus) {
  // Normal constructor with name and description
  frAdcPedRaw       = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseIntRaw  = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseAmpRaw  = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseTimeRaw = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPed          = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseInt     = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseAmp     = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  frAdcPulseTime    = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);
  fAdcErrorFlag     = new TClonesArray("THcSignalHit", MaxNumCerPmt * MaxNumAdcPulse);

  fNumAdcHits         = vector<Int_t>(MaxNumCerPmt, 0.0);
  fNumGoodAdcHits     = vector<Int_t>(MaxNumCerPmt, 0.0);
  fNumTracksMatched   = vector<Int_t>(MaxNumCerPmt, 0.0);
  fNumTracksFired     = vector<Int_t>(MaxNumCerPmt, 0.0);
  fNpe                = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcPed         = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcMult        = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcHitUsed     = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcPulseInt    = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcPulseIntRaw = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcPulseAmp    = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcPulseTime   = vector<Double_t>(MaxNumCerPmt, 0.0);
  fGoodAdcTdcDiffTime = vector<Double_t>(MaxNumCerPmt, 0.0);

  InitArrays();
}

//_____________________________________________________________________________
THcCherenkov::THcCherenkov() {
  // Constructor
  frAdcPedRaw       = NULL;
  frAdcPulseIntRaw  = NULL;
  frAdcPulseAmpRaw  = NULL;
  frAdcPulseTimeRaw = NULL;
  frAdcPed          = NULL;
  frAdcPulseInt     = NULL;
  frAdcPulseAmp     = NULL;
  frAdcPulseTime    = NULL;
  fAdcErrorFlag     = NULL;

  InitArrays();
}

//_____________________________________________________________________________
THcCherenkov::~THcCherenkov() {
  // Destructor
  delete frAdcPedRaw;
  frAdcPedRaw = NULL;
  delete frAdcPulseIntRaw;
  frAdcPulseIntRaw = NULL;
  delete frAdcPulseAmpRaw;
  frAdcPulseAmpRaw = NULL;
  delete frAdcPulseTimeRaw;
  frAdcPulseTimeRaw = NULL;
  delete frAdcPed;
  frAdcPed = NULL;
  delete frAdcPulseInt;
  frAdcPulseInt = NULL;
  delete frAdcPulseAmp;
  frAdcPulseAmp = NULL;
  delete frAdcPulseTime;
  frAdcPulseTime = NULL;
  delete fAdcErrorFlag;
  fAdcErrorFlag = NULL;

  DeleteArrays();
}

//_____________________________________________________________________________
void THcCherenkov::InitArrays() {
  fGain     = NULL;
  fPedSum   = NULL;
  fPedSum2  = NULL;
  fPedLimit = NULL;
  fPedMean  = NULL;
  fPedCount = NULL;
  fPed      = NULL;
  fThresh   = NULL;
}
//_____________________________________________________________________________
void THcCherenkov::DeleteArrays() {
  delete[] fGain;
  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;
}

//_____________________________________________________________________________
THaAnalysisObject::EStatus THcCherenkov::Init(const TDatime& date) {
  // cout << "THcCherenkov::Init for: " << GetName() << endl;

  string EngineDID = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
  std::transform(EngineDID.begin(), EngineDID.end(), EngineDID.begin(), ::toupper);
  if (gHcDetectorMap->FillMap(fDetMap, EngineDID.c_str()) < 0) {
    static const char* const here = "Init()";
    Error(Here(here), "Error filling detectormap for %s.", EngineDID.c_str());
    return kInitError;
  }

  EStatus status;
  if ((status = THaNonTrackingDetector::Init(date)))
    return fStatus = status;

  // Should probably put this in ReadDatabase as we will know the
  // maximum number of hits after setting up the detector map
  InitHitList(fDetMap, "THcCherenkovHit", fDetMap->GetTotNumChan() + 1, 0, fADC_RefTimeCut);

  THcHallCSpectrometer* app = dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
  if (!app || !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector("hod")))) {
    static const char* const here = "ReadDatabase()";
    Warning(Here(here), "Hodoscope \"%s\" not found. ", "hod");
  }

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

//_____________________________________________________________________________
Int_t THcCherenkov::ReadDatabase(const TDatime& date) {
  // This function is called by THaDetectorBase::Init() once at the beginning
  // of the analysis.

  // cout << "THcCherenkov::ReadDatabase for: " << GetName() << endl; // Ahmed

  string prefix = string(GetApparatus()->GetName()).substr(0, 1) + GetName();
  std::transform(prefix.begin(), prefix.end(), prefix.begin(), ::tolower);

  CreateMissReportParms(prefix.c_str());

  fNRegions = 4; // Defualt if not set in paramter file

  DBRequest list_1[] = {{"_tot_pmts", &fNelem, kInt}, {"_num_regions", &fNRegions, kInt}, {0}};

  gHcParms->LoadParmValues(list_1, prefix.c_str());

  Bool_t optional = true;

  _det_logger->info("Created Cherenkov detector {}.{} with {} PMTs", GetApparatus()->GetName(),
                    GetName(), fNelem);
  // cout << "Created Cherenkov detector " << GetApparatus()->GetName() << "."
  //     << 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];
  // Region parameters
  fRegionsValueMax = fNRegions * 8;
  fRegionValue     = new Double_t[fRegionsValueMax];

  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.;
  }
  fDebugAdc       = 0; // Set ADC debug parameter to false unless set in parameter file
  fAdcTdcOffset   = 0.0;
  fADC_RefTimeCut = 0;

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

  // if (fDebugAdc) cout << "Cherenkov ADC Debug Flag Set To TRUE" << endl;

  fIsInit = true;

  // cout << "Track Matching Parameters for: " << GetName() << endl;
  // for (Int_t iregion = 0; iregion < fNRegions; iregion++) {
  //   cout << "Region = " << iregion + 1 << endl;
  //   for (Int_t ivalue = 0; ivalue < 8; ivalue++)
  //     cout << fRegionValue[GetIndex(iregion, ivalue)] << "  ";
  //   cout << endl;
  // }

  // Create arrays to hold pedestal results
  InitializePedestals();

  return kOK;
}

//_____________________________________________________________________________
Int_t THcCherenkov::DefineVariables(EMode mode) {
  // Initialize global variables for histogramming and tree
  // cout << "THcCherenkov::DefineVariables called for: " << GetName() << endl;

  if (mode == kDefine && fIsSetup)
    return kOK;
  fIsSetup = (mode == kDefine);

  // Register variables in global list

  if (fDebugAdc) {
    RVarDef vars[] = {
        {"numAdcHits", "Number of ADC Hits Per PMT", "fNumAdcHits"},     // Cherenkov occupancy
        {"totNumAdcHits", "Total Number of ADC Hits", "fTotNumAdcHits"}, // Cherenkov multiplicity
        {"adcPedRaw", "Raw ADC pedestals", "frAdcPedRaw.THcSignalHit.GetData()"},
        {"adcPulseIntRaw", "Raw ADC pulse integrals", "frAdcPulseIntRaw.THcSignalHit.GetData()"},
        {"adcPulseAmpRaw", "Raw ADC pulse amplitudes", "frAdcPulseAmpRaw.THcSignalHit.GetData()"},
        {"adcPulseTimeRaw", "Raw ADC pulse times", "frAdcPulseTimeRaw.THcSignalHit.GetData()"},
        {"adcPed", "ADC pedestals", "frAdcPed.THcSignalHit.GetData()"},
        {"adcPulseInt", "ADC pulse integrals", "frAdcPulseInt.THcSignalHit.GetData()"},
        {"adcPulseAmp", "ADC pulse amplitudes", "frAdcPulseAmp.THcSignalHit.GetData()"},
        {"adcPulseTime", "ADC pulse times", "frAdcPulseTime.THcSignalHit.GetData()"},
        {0}};
    DefineVarsFromList(vars, 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}};

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

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

void THcCherenkov::Clear(Option_t* opt) {
  // Clear the hit lists
  fNhits               = 0;
  fTotNumAdcHits       = 0;
  fTotNumGoodAdcHits   = 0;
  fTotNumTracksMatched = 0;
  fTotNumTracksFired   = 0;

  fXAtCer = 0.0;
  fYAtCer = 0.0;

  fNpeSum = 0.0;

  frAdcPedRaw->Clear();
  frAdcPulseIntRaw->Clear();
  frAdcPulseAmpRaw->Clear();
  frAdcPulseTimeRaw->Clear();

  frAdcPed->Clear();
  frAdcPulseInt->Clear();
  frAdcPulseAmp->Clear();
  frAdcPulseTime->Clear();
  fAdcErrorFlag->Clear();

  for (UInt_t ielem = 0; ielem < fNumAdcHits.size(); ielem++)
    fNumAdcHits.at(ielem) = 0;
  for (UInt_t ielem = 0; ielem < fNumGoodAdcHits.size(); ielem++)
    fNumGoodAdcHits.at(ielem) = 0;
  for (UInt_t ielem = 0; ielem < fNumTracksMatched.size(); ielem++)
    fNumTracksMatched.at(ielem) = 0;
  for (UInt_t ielem = 0; ielem < fNumTracksFired.size(); ielem++)
    fNumTracksFired.at(ielem) = 0;
  for (UInt_t ielem = 0; ielem < fGoodAdcPed.size(); ielem++) {
    fGoodAdcPed.at(ielem)         = 0.0;
    fGoodAdcMult.at(ielem)        = 0.0;
    fGoodAdcHitUsed.at(ielem)     = 0.0;
    fGoodAdcPulseInt.at(ielem)    = 0.0;
    fGoodAdcPulseIntRaw.at(ielem) = 0.0;
    fGoodAdcPulseAmp.at(ielem)    = 0.0;
    fGoodAdcPulseTime.at(ielem)   = kBig;
    fGoodAdcTdcDiffTime.at(ielem) = kBig;
    fNpe.at(ielem)                = 0.0;
  }
}

//_____________________________________________________________________________

Int_t THcCherenkov::Decode(const THaEvData& evdata) {
  // 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);

  if (gHaCuts->Result("Pedestal_event")) {
    AccumulatePedestals(fRawHitList);
    fAnalyzePedestals = 1; // Analyze pedestals first normal events
    return (0);
  }

  if (fAnalyzePedestals) {
    CalculatePedestals();
    fAnalyzePedestals = 0; // Don't analyze pedestals next event
  }

  Int_t  ihit      = 0;
  UInt_t nrAdcHits = 0;

  while (ihit < fNhits) {

    THcCherenkovHit* hit       = (THcCherenkovHit*)fRawHitList->At(ihit);
    Int_t            npmt      = hit->fCounter;
    THcRawAdcHit&    rawAdcHit = hit->GetRawAdcHitPos();

    for (UInt_t thit = 0; thit < rawAdcHit.GetNPulses(); thit++) {

      ((THcSignalHit*)frAdcPedRaw->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPedRaw());
      ((THcSignalHit*)frAdcPed->ConstructedAt(nrAdcHits))->Set(npmt, rawAdcHit.GetPed());

      ((THcSignalHit*)frAdcPulseIntRaw->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseIntRaw(thit));
      ((THcSignalHit*)frAdcPulseInt->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseInt(thit));

      ((THcSignalHit*)frAdcPulseAmpRaw->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseAmpRaw(thit));
      ((THcSignalHit*)frAdcPulseAmp->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseAmp(thit));

      ((THcSignalHit*)frAdcPulseTimeRaw->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseTimeRaw(thit));
      ((THcSignalHit*)frAdcPulseTime->ConstructedAt(nrAdcHits))
          ->Set(npmt, rawAdcHit.GetPulseTime(thit) + fAdcTdcOffset);

      if (rawAdcHit.GetPulseAmpRaw(thit) > 0)
        ((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 0);
      if (rawAdcHit.GetPulseAmpRaw(thit) <= 0)
        ((THcSignalHit*)fAdcErrorFlag->ConstructedAt(nrAdcHits))->Set(npmt, 1);

      ++nrAdcHits;
      fTotNumAdcHits++;
      fNumAdcHits.at(npmt - 1) = npmt;
    }
    ihit++;
  }
  return ihit;
}

//_____________________________________________________________________________
Int_t THcCherenkov::ApplyCorrections(void) { return (0); }

//_____________________________________________________________________________
Int_t THcCherenkov::CoarseProcess(TClonesArray&) {
  Double_t StartTime = 0.0;
  if (fglHod)
    StartTime = fglHod->GetStartTime();

  // Loop over the elements in the TClonesArray
  for (Int_t ielem = 0; ielem < frAdcPulseInt->GetEntries(); ielem++) {

    Int_t    npmt     = ((THcSignalHit*)frAdcPulseInt->ConstructedAt(ielem))->GetPaddleNumber() - 1;
    Double_t pulsePed = ((THcSignalHit*)frAdcPed->ConstructedAt(ielem))->GetData();
    Double_t pulseInt = ((THcSignalHit*)frAdcPulseInt->ConstructedAt(ielem))->GetData();
    Double_t pulseIntRaw    = ((THcSignalHit*)frAdcPulseIntRaw->ConstructedAt(ielem))->GetData();
    Double_t pulseAmp       = ((THcSignalHit*)frAdcPulseAmp->ConstructedAt(ielem))->GetData();
    Double_t pulseTime      = ((THcSignalHit*)frAdcPulseTime->ConstructedAt(ielem))->GetData();
    Double_t adctdcdiffTime = StartTime - pulseTime;
    Bool_t   errorFlag      = ((THcSignalHit*)fAdcErrorFlag->ConstructedAt(ielem))->GetData();
    Bool_t   pulseTimeCut =
        adctdcdiffTime > fAdcTimeWindowMin[npmt] && adctdcdiffTime < fAdcTimeWindowMax[npmt];

    if (!errorFlag) {
      fGoodAdcMult.at(npmt) += 1;
    }

    // By default, the last hit within the timing cut will be considered "good"
    if (!errorFlag && pulseTimeCut) {
      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);
      fNpeSum += fNpe.at(npmt);

      fTotNumGoodAdcHits++;
      fNumGoodAdcHits.at(npmt) = npmt + 1;
    }
  }
  return 0;
}

//_____________________________________________________________________________
Int_t THcCherenkov::FineProcess(TClonesArray& tracks) {

  Int_t nTracks = tracks.GetLast() + 1;

  for (Int_t itrack = 0; itrack < nTracks; itrack++) {

    THaTrack* track = dynamic_cast<THaTrack*>(tracks[itrack]);
    if (track->GetIndex() != 0)
      continue; // Select the best track

    Double_t trackChi2    = track->GetChi2();
    Int_t    trackNDoF    = track->GetNDoF();
    Double_t trackRedChi2 = trackChi2 / trackNDoF;
    Double_t trackBeta    = track->GetBeta();
    Double_t trackEnergy  = track->GetEnergy();
    Double_t trackMom     = track->GetP();
    Double_t trackDp      = track->GetDp();
    Double_t trackENorm   = trackEnergy / trackMom;
    Double_t trackXfp     = track->GetX();
    Double_t trackYfp     = track->GetY();
    Double_t trackTheta   = track->GetTheta();
    Double_t trackPhi     = track->GetPhi();

    Bool_t trackRedChi2Cut = trackRedChi2 > fRedChi2Min && trackRedChi2 < fRedChi2Max;
    Bool_t trackBetaCut    = trackBeta > fBetaMin && trackBeta < fBetaMax;
    Bool_t trackENormCut   = trackENorm > fENormMin && trackENorm < fENormMax;
    Bool_t trackDpCut      = trackDp > fDpMin && trackDp < fDpMax;
    fXAtCer                = trackXfp + trackTheta * fMirrorZPos;
    fYAtCer                = trackYfp + trackPhi * fMirrorZPos;

    if (trackRedChi2Cut && trackBetaCut && trackENormCut && trackDpCut) {

      // Project the track to the Cherenkov mirror planes

      // cout << "Cherenkov Detector: " << GetName() << " has fNRegions = " << fNRegions << endl;
      // cout << "nTracks = " << nTracks << "\t" << "trackChi2 = " << trackChi2
      // 	   << "\t" << "trackNDof = " << trackNDoF << "\t" << "trackRedChi2 = " <<
      // trackRedChi2 << endl; cout << "trackBeta = " << trackBeta << "\t" << "trackEnergy = " <<
      // trackEnergy << "\t"
      // 	   << "trackMom = " << trackMom << "\t" << "trackENorm = " << trackENorm << endl;
      // cout << "trackXfp = " << trackXfp << "\t" << "trackYfp = " << trackYfp << "\t"
      // 	   << "trackTheta = " << trackTheta << "\t" << "trackPhi = " << trackPhi << endl;
      // cout << "fMirrorZPos = " << fMirrorZPos << "\t" << "fXAtCer = " << fXAtCer << "\t" <<
      // "fYAtCer = " << fYAtCer << endl; cout <<
      // "=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:=:" << endl;

      for (Int_t iregion = 0; iregion < fNRegions; iregion++) {

        if ((TMath::Abs(fRegionValue[GetIndex(iregion, 0)] - fXAtCer) <
             fRegionValue[GetIndex(iregion, 4)]) &&
            (TMath::Abs(fRegionValue[GetIndex(iregion, 1)] - fYAtCer) <
             fRegionValue[GetIndex(iregion, 5)]) &&
            (TMath::Abs(fRegionValue[GetIndex(iregion, 2)] - trackTheta) <
             fRegionValue[GetIndex(iregion, 6)]) &&
            (TMath::Abs(fRegionValue[GetIndex(iregion, 3)] - trackPhi) <
             fRegionValue[GetIndex(iregion, 7)])) {

          fTotNumTracksMatched++;
          fNumTracksMatched.at(iregion) = iregion + 1;

          if (fNpeSum > fNpeThresh) {
            fTotNumTracksFired++;
            fNumTracksFired.at(iregion) = iregion + 1;
          } // NPE threshold cut
        }   // Regional cuts
      }     // Loop over regions
    }       // Tracking cuts
  }         // Track loop

  return 0;
}

//_____________________________________________________________________________
void THcCherenkov::InitializePedestals() {
  fNPedestalEvents = 0;
  fMinPeds         = 500; // In engine, this is set in parameter file
  fPedSum          = new Int_t[fNelem];
  fPedSum2         = new Int_t[fNelem];
  fPedCount        = new Int_t[fNelem];
  fPed             = new Double_t[fNelem];
  fThresh          = new Double_t[fNelem];
  for (Int_t i = 0; i < fNelem; i++) {
    fPedSum[i]   = 0;
    fPedSum2[i]  = 0;
    fPedCount[i] = 0;
  }
}

//_____________________________________________________________________________
void THcCherenkov::AccumulatePedestals(TClonesArray* rawhits) {
  // Extract data from the hit list, accumulating into arrays for
  // calculating pedestals

  Int_t nrawhits = rawhits->GetLast() + 1;

  Int_t ihit = 0;
  while (ihit < nrawhits) {
    THcCherenkovHit* hit = (THcCherenkovHit*)rawhits->At(ihit);

    Int_t element = hit->fCounter - 1;
    Int_t nadc    = hit->GetRawAdcHitPos().GetPulseIntRaw();
    if (nadc <= fPedLimit[element]) {
      fPedSum[element] += nadc;
      fPedSum2[element] += nadc * nadc;
      fPedCount[element]++;
      if (fPedCount[element] == fMinPeds / 5) {
        fPedLimit[element] = 100 + fPedSum[element] / fPedCount[element];
      }
    }
    ihit++;
  }

  fNPedestalEvents++;

  return;
}

//_____________________________________________________________________________
void THcCherenkov::CalculatePedestals() {
  // Use the accumulated pedestal data to calculate pedestals
  // Later add check to see if pedestals have drifted ("Danger Will Robinson!")
  //  cout << "Plane: " << fPlaneNum << endl;
  for (Int_t i = 0; i < fNelem; i++) {

    // PMT tubes
    fPed[i]    = ((Double_t)fPedSum[i]) / TMath::Max(1, fPedCount[i]);
    fThresh[i] = fPed[i] + 15;

    // Just a copy for now, but allow the possibility that fXXXPedMean is set
    // in a parameter file and only overwritten if there is a sufficient number of
    // pedestal events.  (So that pedestals are sensible even if the pedestal events were
    // not acquired.)
    if (fMinPeds > 0) {
      if (fPedCount[i] > fMinPeds) {
        fPedMean[i] = fPed[i];
      }
    }
  }
  //  cout << " " << endl;
}

//_____________________________________________________________________________
Int_t THcCherenkov::GetIndex(Int_t nRegion, Int_t nValue) { return fNRegions * nValue + nRegion; }

//_____________________________________________________________________________
void THcCherenkov::Print(const Option_t* opt) const {
  THaNonTrackingDetector::Print(opt);
  // Print out the pedestals
  cout << endl;
  cout << "Cherenkov Pedestals" << endl;
  // Ahmed
  cout << "No.   ADC" << endl;
  for (Int_t i = 0; i < fNelem; i++)
    cout << " " << i << "    " << fPed[i] << endl;
  cout << endl;
}

//_____________________________________________________________________________
Double_t THcCherenkov::GetCerNPE() { return fNpeSum; }
//_____________________________________________________________________________
Int_t THcCherenkov::End(THaRunBase* run) {
  MissReport(Form("%s.%s", GetApparatus()->GetName(), GetName()));
  return 0;
}
ClassImp(THcCherenkov)
    ////////////////////////////////////////////////////////////////////////////////