From 39317b13faf07742e66a5c7156205bcae4d83767 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Tue, 6 Feb 2018 12:08:28 -0500
Subject: [PATCH] Put a minimum cut on reference times.    First hit above this
 cut is taken as the reference time.    Cuts are per detector and read as an
 optional parameter in    each detector class.  TDCs and Flash ADC times have
 separate cuts.

Implement a second method of selecting best reference time.
  If InitHitList is passed a possitive reference time cut,
    use the first reference time above that cut.  If none of the
    hits pass the cut, then no reference time is set
  If InitHitList is passed a negative cut, the cut is taken to be
    the absolute value of the cut, but a reference time is garuanteed
    to be set as long as there is at least one reference time hit.
    If none of the reference time hits pass the cut, then the last
    hit is taken to be the reference time.  (Should be the largest.)

Allow reference time cuts for all the spectrometer detector classes
  Hodoscopes, drift chambers, Aerogels, Gas Cherenkovs, Calorimeters
---
 src/THcAerogel.cxx   |  14 +++--
 src/THcAerogel.h     |   2 +
 src/THcCherenkov.cxx |  11 ++--
 src/THcCherenkov.h   |   2 +
 src/THcDC.cxx        |   7 ++-
 src/THcDC.h          |   2 +
 src/THcHitList.cxx   | 142 +++++++++++++++++++++++++++++--------------
 src/THcHitList.h     |   7 ++-
 src/THcHodoscope.cxx |   7 ++-
 src/THcHodoscope.h   |   3 +
 src/THcShower.cxx    |   5 +-
 11 files changed, 140 insertions(+), 62 deletions(-)

diff --git a/src/THcAerogel.cxx b/src/THcAerogel.cxx
index 9f9e421..d43d30b 100644
--- a/src/THcAerogel.cxx
+++ b/src/THcAerogel.cxx
@@ -185,14 +185,15 @@ THaAnalysisObject::EStatus THcAerogel::Init( const TDatime& date )
     return kInitError;
   }
 
-  // Should probably put this in ReadDatabase as we will know the
-  // maximum number of hits after setting up the detector map
-  InitHitList(fDetMap, "THcAerogelHit", fDetMap->GetTotNumChan()+1);
-
   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, "THcAerogelHit", fDetMap->GetTotNumChan()+1,
+	      0, fADC_RefTimeCut);
+
  THcHallCSpectrometer *app=dynamic_cast<THcHallCSpectrometer*>(GetApparatus());
    if(  !app ||
       !(fglHod = dynamic_cast<THcHodoscope*>(app->GetDetector("hod"))) ) {
@@ -333,13 +334,14 @@ Int_t THcAerogel::ReadDatabase( const TDatime& date )
     {"aero_tdc_offset",       &fTdcOffset,        kInt,    0,               optional},
     {"aero_min_peds",         &fMinPeds,          kInt,    0,               optional},
     {"aero_region",           &fRegionValue[0],   kDouble, (UInt_t) fRegionsValueMax},
-
+    {"aero_adcrefcut",        &fADC_RefTimeCut,   kInt,    0, 1},
     {0}
   };
 
   fSixGevData = 0; // Set 6 GeV data parameter to false unless set in parameter file
   fDebugAdc   = 0; // Set ADC debug parameter to false unless set in parameter file
   fAdcTdcOffset = 0.0;
+  fADC_RefTimeCut = 0;
 
   gHcParms->LoadParmValues((DBRequest*)&list, prefix);
 
diff --git a/src/THcAerogel.h b/src/THcAerogel.h
index 34be99e..10f1a1d 100644
--- a/src/THcAerogel.h
+++ b/src/THcAerogel.h
@@ -44,6 +44,8 @@ class THcAerogel : public THaNonTrackingDetector, public THcHitList {
   Int_t fNhits;
   Bool_t* fPresentP;
 
+  Int_t fADC_RefTimeCut;
+
   // 12 GeV variables
   // Vector/TClonesArray length parameters
   static const Int_t MaxNumAdcPulse   = 4;
diff --git a/src/THcCherenkov.cxx b/src/THcCherenkov.cxx
index 1e44e79..82eed26 100644
--- a/src/THcCherenkov.cxx
+++ b/src/THcCherenkov.cxx
@@ -148,14 +148,15 @@ THaAnalysisObject::EStatus THcCherenkov::Init( const TDatime& date )
     return kInitError;
   }
 
-  // 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);
-
   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"))) ) {
@@ -224,11 +225,13 @@ Int_t THcCherenkov::ReadDatabase( const TDatime& date )
     {"_adcTimeWindowMax", &fAdcTimeWindowMax, kDouble},
     {"_adc_tdc_offset",   &fAdcTdcOffset,     kDouble, 0, 1},
     {"_region",           &fRegionValue[0],   kDouble,  (UInt_t) fRegionsValueMax},
+    {"_adcrefcut",        &fADC_RefTimeCut,   kInt,    0, 1},
     {0}
   };
 
   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());
 
diff --git a/src/THcCherenkov.h b/src/THcCherenkov.h
index 5c9e0cc..6b852ec 100644
--- a/src/THcCherenkov.h
+++ b/src/THcCherenkov.h
@@ -51,6 +51,8 @@ class THcCherenkov : public THaNonTrackingDetector, public THcHitList {
   Int_t     fDebugAdc;
   Double_t* fWidth;
 
+  Int_t     fADC_RefTimeCut;
+
   Int_t     fNhits;
   Int_t     fTotNumAdcHits;
   Int_t     fTotNumGoodAdcHits;
diff --git a/src/THcDC.cxx b/src/THcDC.cxx
index 44ab6bc..b4b7d91 100644
--- a/src/THcDC.cxx
+++ b/src/THcDC.cxx
@@ -117,10 +117,12 @@ void THcDC::Setup(const char* name, const char* description)
     {"dc_tdc_time_per_channel",&fNSperChan, kDouble},
     {"dc_wire_velocity",&fWireVelocity,kDouble},
     {"dc_plane_names",&planenamelist, kString},
-	{"dc_version", &fVersion, kInt, 0, optional},
+    {"dc_version", &fVersion, kInt, 0, optional},
+    {"dc_tdcrefcut", &fTDC_RefTimeCut, kInt, 0, 1},
     {0}
   };
 
+  fTDC_RefTimeCut = 0;		// Minimum allowed reference times
   gHcParms->LoadParmValues((DBRequest*)&list,fPrefix);
 
   if(fVersion==0) {
@@ -207,7 +209,8 @@ THaAnalysisObject::EStatus THcDC::Init( const TDatime& date )
 
   // Should probably put this in ReadDatabase as we will know the
   // maximum number of hits after setting up the detector map
-  InitHitList(fDetMap, "THcRawDCHit", fDetMap->GetTotNumChan()+1);
+  InitHitList(fDetMap, "THcRawDCHit", fDetMap->GetTotNumChan()+1,
+	      fTDC_RefTimeCut, 0);
 
   EStatus status;
   // This triggers call of ReadDatabase and DefineVariables
diff --git a/src/THcDC.h b/src/THcDC.h
index 7d1c6ae..43667b8 100644
--- a/src/THcDC.h
+++ b/src/THcDC.h
@@ -93,6 +93,8 @@ protected:
   Int_t fdebugprintdecodeddc;
   Int_t fHMSStyleChambers;
 
+  Int_t fTDC_RefTimeCut;
+
   UInt_t fNDCTracks;
   TClonesArray* fDCTracks;     // Tracks found from stubs (THcDCTrack obj)
   // Calibration
diff --git a/src/THcHitList.cxx b/src/THcHitList.cxx
index e934490..3489065 100644
--- a/src/THcHitList.cxx
+++ b/src/THcHitList.cxx
@@ -42,13 +42,29 @@ initialize a hit array of hits of class hitclass
 */
 
 void THcHitList::InitHitList(THaDetMap* detmap,
-				  const char *hitclass, Int_t maxhits) {
-
-
+			     const char *hitclass, Int_t maxhits,
+			     Int_t tdcref_cut, Int_t adcref_cut) {
+  cout << "InitHitList: " << hitclass << " RefTimeCuts: " << tdcref_cut << " " << adcref_cut << endl;
   fRawHitList = new TClonesArray(hitclass, maxhits);
   fRawHitClass = fRawHitList->GetClass();
   fNMaxRawHits = maxhits;
   fNRawHits = 0;
+
+  if(tdcref_cut >= 0) {
+    fTDC_RefTimeBest = kFALSE;
+    fTDC_RefTimeCut = tdcref_cut;
+  } else {
+    fTDC_RefTimeBest = kTRUE;
+    fTDC_RefTimeCut = -tdcref_cut;
+  }
+  if(adcref_cut >= 0) {
+    fADC_RefTimeBest = kFALSE;
+    fADC_RefTimeCut = adcref_cut;
+  } else {
+    fADC_RefTimeBest = kTRUE;
+    fADC_RefTimeCut = -adcref_cut;
+  }
+
   for(Int_t i=0;i<maxhits;i++) {
     fRawHitList->ConstructedAt(i);
   }
@@ -163,36 +179,52 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
       if(evdata.IsMultifunction(fRefIndexMaps[i].crate,
 				fRefIndexMaps[i].slot)) { // Multifunction module (e.g. FADC)
 	// Make sure at least one pulse
-	if(evdata.GetNumEvents(Decoder::kPulseTime, fRefIndexMaps[i].crate,
-			       fRefIndexMaps[i].slot,
-			       fRefIndexMaps[i].channel) > 0) {
-	  fRefIndexMaps[i].reftime =
-	    evdata.GetData(Decoder::kPulseTime,fRefIndexMaps[i].crate,fRefIndexMaps[i].slot,
-			   fRefIndexMaps[i].channel,0);
-	  fRefIndexMaps[i].hashit = kTRUE;
-	} else {
-	  fRefIndexMaps[i].hashit = kFALSE;
+	Int_t nrefhits = evdata.GetNumEvents(Decoder::kPulseTime,
+					     fRefIndexMaps[i].crate,
+					     fRefIndexMaps[i].slot,
+					     fRefIndexMaps[i].channel);
+	fRefIndexMaps[i].hashit = kFALSE;
+	Bool_t goodreftime=kFALSE;
+	Int_t reftime = 0;
+	for(Int_t ihit=0; ihit<nrefhits; ihit++) {
+	  reftime = evdata.GetData(Decoder::kPulseTime,fRefIndexMaps[i].crate,
+				   fRefIndexMaps[i].slot, fRefIndexMaps[i].channel,ihit);
+	  if(reftime >= fADC_RefTimeCut) {
+	    goodreftime = kTRUE;
+	    break;
+	  }
 	}
-      } else {
-	if(evdata.GetNumHits(fRefIndexMaps[i].crate,
-			     fRefIndexMaps[i].slot,
-			     fRefIndexMaps[i].channel) > 0) {
-	  // Only take first hit in this reference channel
-	  // Here we need to check if it is a multifunction type and get the time
-	  // if it is.
-	  fRefIndexMaps[i].reftime =
-	    evdata.GetData(fRefIndexMaps[i].crate,fRefIndexMaps[i].slot,
-			   fRefIndexMaps[i].channel,0);
+	if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
+	  fRefIndexMaps[i].reftime = reftime;
 	  fRefIndexMaps[i].hashit = kTRUE;
-	} else {
-	  fRefIndexMaps[i].hashit = kFALSE;
+	}
+      } else {			// Assume this is a TDC
+	Int_t nrefhits = evdata.GetNumHits(fRefIndexMaps[i].crate,
+					   fRefIndexMaps[i].slot,
+					   fRefIndexMaps[i].channel);
+	fRefIndexMaps[i].hashit = kFALSE;
+	// Only take first hit in this reference channel that is bigger
+	// then fTDC_RefTimeCut
+	Bool_t goodreftime=kFALSE;
+	Int_t reftime = 0;
+	for(Int_t ihit=0; ihit<nrefhits; ihit++) {
+	  reftime = evdata.GetData(fRefIndexMaps[i].crate,fRefIndexMaps[i].slot,
+				   fRefIndexMaps[i].channel,ihit);
+	  if(reftime >= fTDC_RefTimeCut) {
+	    goodreftime = kTRUE;
+	    break;
+	  }
+	}
+	if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
+	    fRefIndexMaps[i].reftime = reftime;
+	    fRefIndexMaps[i].hashit = kTRUE;
 	}
       }
     }
   }
   for ( Int_t i=0; i < fdMap->GetSize(); i++ ) {
     THaDetMap::Module* d = fdMap->GetModule(i);
-
+    
     // Loop over all channels that have a hit.
     //    cout << "Crate/Slot: " << d->crate << "/" << d->slot << endl;
     Int_t plane = d->plane;
@@ -243,20 +275,27 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	  // cout << "Signal " << signal << "=" << data << endl;
 	  rawhit->SetData(signal,data);
 	}
-	// Get the reference time.  Only take the first hit
-	// If a reference channel
-	// was specified, it takes precidence of reference index
+	// Get the reference time.
 	if(d->refchan >= 0) {
-	  if( evdata.GetNumHits(d->crate,d->slot,d->refchan) > 0) {
-	    Int_t reftime = evdata.GetData(d->crate, d->slot, d->refchan, 0);
+	  Int_t nrefhits = evdata.GetNumHits(d->crate,d->slot,d->refchan);
+	  Bool_t goodreftime=kFALSE;
+	  Int_t reftime=0;
+	  for(Int_t ihit=0; ihit<nrefhits; ihit++) {
+	    reftime = evdata.GetData(d->crate, d->slot, d->refchan, ihit);
+	    if(reftime >= fTDC_RefTimeCut) {
+	      goodreftime = kTRUE;
+	      break;
+	    }
+	  }
+	  // If RefTimeBest flag set, take the last hit if none of the
+	  // hits make the RefTimeCut
+	  if(goodreftime || (nrefhits>0 && fTDC_RefTimeBest)) {
 	    rawhit->SetReference(signal, reftime);
-	  } else {
-	    if(!suppresswarnings) {
-	      cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
-		" missing for (" << d->crate << ", " << d->slot <<
-		", " << chan << ")" << endl;
+	  } else if (!suppresswarnings) {
+	    cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
+	      " missing for (" << d->crate << ", " << d->slot <<
+	      ", " << chan << ")" << endl;
 	      tdcref_miss = kTRUE;
-	    }
 	  }
 	} else {
 	  if(d->refindex >=0 && d->refindex < fNRefIndex) {
@@ -284,10 +323,10 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	    fNPED = fPSE125->GetNPED(d->crate);
 	    fHaveFADCInfo = kTRUE;
 	  }
-// Set F250 parameters.
+	  // Set F250 parameters.
           rawhit->SetF250Params(fNSA, fNSB, fNPED);
         }
-
+	
 	// Copy the samples
 	Int_t nsamples=evdata.GetNumEvents(Decoder::kSampleADC, d->crate, d->slot, chan);
 
@@ -309,16 +348,26 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 	}
 	// Get the reference time for the FADC pulse time
 	if(d->refchan >= 0) {	// Reference time for the slot
-	  if(evdata.GetNumEvents(Decoder::kPulseIntegral, d->crate, d->slot, d->refchan) > 0) {
-	    Int_t reftime = evdata.GetData(Decoder::kPulseTime, d->crate, d->slot, d->refchan, 0);
+	  Int_t nrefhits = evdata.GetNumEvents(Decoder::kPulseIntegral,
+					       d->crate, d->slot, d->refchan);
+	  Bool_t goodreftime=kFALSE;
+	  Int_t reftime = 0;
+	  for(Int_t ihit=0; ihit<nrefhits; ihit++) {
+	    reftime = evdata.GetData(Decoder::kPulseTime, d->crate, d->slot, d->refchan, ihit);
+	    if(reftime >= fADC_RefTimeCut) {
+	      goodreftime=kTRUE;
+	      break;
+	    }
+	  }
+	  // If RefTimeBest flag set, take the last hit if none of the
+	  // hits make the RefTimeCut
+	  if(goodreftime || (nrefhits>0 && fADC_RefTimeBest)) {
 	    rawhit->SetReference(signal, reftime);
-	  } else {
-	    if(!suppresswarnings) {
-	      cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
-		" missing for (" << d->crate << ", " << d->slot <<
-		", " << chan << ")" << endl;
+	  } else if (!suppresswarnings) {
+	    cout << "HitList(event=" << evdata.GetEvNum() << "): refchan " << d->refchan <<
+	      " missing for (" << d->crate << ", " << d->slot <<
+	      ", " << chan << ")" << endl;
 	      adcref_miss = kTRUE;
-	    }
 	  }
 	} else {
 	  if(d->refindex >=0 && d->refindex < fNRefIndex) {
@@ -344,7 +393,6 @@ Int_t THcHitList::DecodeToHitList( const THaEvData& evdata, Bool_t suppresswarni
 
   fNTDCRef_miss += (tdcref_miss ? 1 : 0);
   fNADCRef_miss += (adcref_miss ? 1 : 0);
-
   return fNRawHits;		// Does anything care what is returned
 }
 void THcHitList::MissReport(const char *name)
diff --git a/src/THcHitList.h b/src/THcHitList.h
index eb2ec85..4b378ce 100644
--- a/src/THcHitList.h
+++ b/src/THcHitList.h
@@ -31,13 +31,18 @@ public:
 
   virtual Int_t DecodeToHitList( const THaEvData&, Bool_t suppress=kFALSE );
   void          InitHitList(THaDetMap* detmap,
-			    const char *hitclass, Int_t maxhits);
+			    const char *hitclass, Int_t maxhits,
+			    Int_t tdcref_cut=0, Int_t adcref_cut=0);
 
   TClonesArray* GetHitList() const {return fRawHitList; }
   void          MissReport(const char *name);
 
   UInt_t         fNRawHits;
   Int_t         fNMaxRawHits;
+  Int_t         fTDC_RefTimeCut;
+  Int_t         fADC_RefTimeCut;
+  Bool_t        fTDC_RefTimeBest;
+  Bool_t        fADC_RefTimeBest;
   TClonesArray* fRawHitList; // List of raw hits
   TClass* fRawHitClass;		  // Class of raw hit object to use
 
diff --git a/src/THcHodoscope.cxx b/src/THcHodoscope.cxx
index fc382e4..3f5e9f9 100644
--- a/src/THcHodoscope.cxx
+++ b/src/THcHodoscope.cxx
@@ -95,9 +95,13 @@ void THcHodoscope::Setup(const char* name, const char* description)
   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);
 
   cout << "Plane Name List : " << planenamelist << endl;
@@ -160,7 +164,8 @@ THaAnalysisObject::EStatus THcHodoscope::Init( const TDatime& date )
   // But it needs to happen before the sub detectors are initialized
   // so that they can get the pointer to the hitlist.
 
-  InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan()+1);
+  InitHitList(fDetMap, "THcRawHodoHit", fDetMap->GetTotNumChan()+1,
+	      fTDC_RefTimeCut, fADC_RefTimeCut);
 
   EStatus status;
   // This triggers call of ReadDatabase and DefineVariables
diff --git a/src/THcHodoscope.h b/src/THcHodoscope.h
index fb9ee76..7edfb57 100644
--- a/src/THcHodoscope.h
+++ b/src/THcHodoscope.h
@@ -122,6 +122,9 @@ protected:
 
   THcCherenkov* fCherenkov;
 
+  Int_t fTDC_RefTimeCut;
+  Int_t fADC_RefTimeCut;
+
   Int_t fAnalyzePedestals;
 
   Int_t fNHits;
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index aefa5fa..a0ebc2e 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -63,9 +63,11 @@ void THcShower::Setup(const char* name, const char* description)
     {"cal_num_layers", &fNLayers, kInt},
     {"cal_layer_names", &layernamelist, kString},
     {"cal_array",&fHasArray, kInt,0, 1},
+    {"cal_adcrefcut", &fADC_RefTimeCut, kInt, 0, 1},
     {0}
   };
 
+  fADC_RefTimeCut = 0;
   gHcParms->LoadParmValues((DBRequest*)&list,prefix);
   fNTotLayers = (fNLayers+(fHasArray!=0?1:0));
   vector<string> layer_names = vsplit(layernamelist);
@@ -141,7 +143,8 @@ THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
   // Should probably put this in ReadDatabase as we will know the
   // maximum number of hits after setting up the detector map
 
-  InitHitList(fDetMap, "THcRawShowerHit", fDetMap->GetTotNumChan()+1);
+  InitHitList(fDetMap, "THcRawShowerHit", fDetMap->GetTotNumChan()+1,
+	      0, fADC_RefTimeCut);
 
   EStatus status;
   if( (status = THaNonTrackingDetector::Init( date )) )
-- 
GitLab