THcCherenkov.cxx

/** \class THcCherenkov
    \ingroup Detectors

\brief Class for gas Cherenkov detectors

*/

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

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


using namespace std;

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

//_____________________________________________________________________________
THcCherenkov::THcCherenkov( const char* name, const char* description,
                            THaApparatus* apparatus ) :
  hcana::ConfigLogging<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( ) :
  hcana::ConfigLogging<THaNonTrackingDetector>()
{
  // 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;
}

//_____________________________________________________________________________
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;

  _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);
}
//_____________________________________________________________________________
inline
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)
////////////////////////////////////////////////////////////////////////////////