From f8f421ce6b42711c3f9574d39e3e57b241181df8 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Wed, 21 Oct 2015 16:37:56 -0400
Subject: [PATCH] Add dynamic pedestal subtraction to THcRawShowerHit class  
 GetData accepts a pedestal range and an integration range.   Samples in the
 integration range are added after subtracting the   average of the pedestal
 range from each sample.   Method also provided to return dynamic pedestal

---
 src/THcRawShowerHit.cxx | 78 ++++++++++++++++++++++++++++++++++++++++-
 src/THcRawShowerHit.h   |  5 ++-
 src/THcShowerArray.cxx  | 22 +++++++++++-
 src/THcShowerArray.h    |  6 ++++
 4 files changed, 108 insertions(+), 3 deletions(-)

diff --git a/src/THcRawShowerHit.cxx b/src/THcRawShowerHit.cxx
index 5c28eac..d509621 100644
--- a/src/THcRawShowerHit.cxx
+++ b/src/THcRawShowerHit.cxx
@@ -6,10 +6,18 @@
 //                                                                           //
 // Contains plane, counter and pos/neg adc and tdc values                    //
 //                                                                           //
+// Enhanced to work with FADC250 data samples.  If fNPosSamples/fNNegSamples //
+// is greater than 1, assume that the data held in the hit is the sampled    //
+// waveform.  Signals 0,1 will return the integrated pulse with dynamic      //
+// pedestal subtraction (first four samples comprise the pedestal).  Signals //
+// 2,3 are reserved for time information.  Signals 4,5 are pedestals and     //
+// 6 and 7 are the straight sum of all the samples.                          //
+//                                                                           //
 ///////////////////////////////////////////////////////////////////////////////
 
 #include "THcRawShowerHit.h"
 #include <iostream>
+#include <cassert>
 
 using namespace std;
 
@@ -22,9 +30,56 @@ void THcRawShowerHit::SetData(Int_t signal, Int_t data) {
   } 
 }
 
+Int_t THcRawShowerHit::GetData(Int_t signal, Int_t isamplow, Int_t isamphigh,
+	      Int_t iintegrallow, Int_t iintegralhigh) {
+  Int_t adcsum=0;
+  Double_t pedestal=0.0;
+
+  if(signal==0 || signal == 1) {
+    pedestal = GetPedestal(signal, isamplow, isamphigh);
+
+#if 0
+    for(Int_t i=0;i<fNPosSamples;i++) {
+      cout << fCounter << " " << i;
+      if(i >= isamplow && i<=isamphigh) {
+	cout << "P ";
+      }	else if (i >= iintegrallow && i<=iintegralhigh) {
+	cout << "D ";
+      } else {
+	cout << "  ";
+      }
+      cout << fADC_pos[i] << " " << fADC_pos[i] - pedestal << endl;
+    }
+#endif
+  }
+  if(signal==4 || signal==5) {
+    pedestal = GetPedestal(signal-4, isamplow, isamphigh);
+    return(pedestal);
+  }
+  if(signal==0 || signal==6) {
+    assert(iintegralhigh<(Int_t) fNPosSamples);
+    for(UInt_t isample=iintegrallow;isample<=iintegralhigh;isample++) {
+      adcsum += fADC_pos[isample] - pedestal;
+    }
+    return(adcsum);
+  } else if (signal==1 || signal==7) {
+    assert(iintegralhigh<(Int_t) fNNegSamples);
+    for(UInt_t isample=iintegrallow;isample<=iintegralhigh;isample++) {
+      adcsum += fADC_neg[isample] - pedestal;
+    }
+    return(adcsum);
+  } 
+  return(-1); // Actually should throw exception
+}
+  
 // Return sum of samples
+// For Fastbus ADC, this will simply be the hardware ADC value as there
+// is just one sample.  For Flash ADCs, this should return the
+// integrated ADC value if the FADC provided that, otherwise the sum
+// of all the samples
 Int_t THcRawShowerHit::GetData(Int_t signal) {
   Int_t adcsum=0;
+
   if(signal==0) {
     for(UInt_t isample=0;isample<fNPosSamples;isample++) {
       adcsum += fADC_pos[isample];
@@ -40,7 +95,7 @@ Int_t THcRawShowerHit::GetData(Int_t signal) {
 }
 
 // Return a requested sample
-Int_t THcRawShowerHit::GetData(Int_t signal, UInt_t isample) {
+Int_t THcRawShowerHit::GetSample(Int_t signal, UInt_t isample) {
   if(signal==0) {
     if(isample >=0 && isample< fNPosSamples) {
       return(fADC_pos[isample]);
@@ -53,6 +108,27 @@ Int_t THcRawShowerHit::GetData(Int_t signal, UInt_t isample) {
   return(-1);
 }
 
+Double_t THcRawShowerHit::GetPedestal(Int_t signal, Int_t isamplow, Int_t isamphigh) {
+  // No error checking on pedestal range
+  Double_t pedestal=0.0;
+  if(signal==0 && fNPosSamples > 1) {
+    if(fNPosSamples > isamphigh) {
+      for(Int_t i = isamplow;i<=isamphigh;i++) {
+	pedestal += fADC_pos[i];
+      }
+      return(pedestal/(isamphigh-isamplow+1));
+    }
+  } else if (signal==1 && fNNegSamples > 1) {
+    if(fNNegSamples > isamphigh) {
+      for(Int_t i = isamplow;i<=isamphigh;i++) {
+	pedestal += fADC_neg[i];
+      }
+      return(pedestal/(isamphigh-isamplow+1));
+    }
+  }
+  return(pedestal);
+}
+      
 // Return the number of samples
 Int_t THcRawShowerHit::GetNSamples(Int_t signal) {
   if(signal==0) {
diff --git a/src/THcRawShowerHit.h b/src/THcRawShowerHit.h
index f7cd000..5a1102d 100644
--- a/src/THcRawShowerHit.h
+++ b/src/THcRawShowerHit.h
@@ -22,7 +22,10 @@ class THcRawShowerHit : public THcRawHit {
 
   void SetData(Int_t signal, Int_t data);
   Int_t GetData(Int_t signal);
-  Int_t GetData(Int_t signal, UInt_t isample);
+  Int_t GetData(Int_t signal, Int_t isamplow, Int_t isamphigh,
+		Int_t iintegrallow, Int_t iintegralhigh);
+  Int_t GetSample(Int_t signal, UInt_t isample);
+  Double_t GetPedestal(Int_t signal, Int_t isamplow, Int_t isamphigh);
   Int_t GetNSamples(Int_t signal);
 
   //  virtual Bool_t  IsSortable () const {return kTRUE; }
diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
index 1d46585..2c57c51 100644
--- a/src/THcShowerArray.cxx
+++ b/src/THcShowerArray.cxx
@@ -46,6 +46,7 @@ THcShowerArray::~THcShowerArray()
   delete fADCHits;
 
   delete [] fA;
+  delete [] fP;
 }
 
 //_____________________________________________________________________________
@@ -77,9 +78,19 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
   prefix[1]='\0';
 
   cout << "Parent name: " << GetParent()->GetPrefix() << endl;
+  fUsingFADC=0;
+  fPedSampLow=0;
+  fPedSampHigh=9;
+  fDataSampLow=23;
+  fDataSampHigh=49;
   DBRequest list[]={
     {"cal_nrows", &fNRows, kInt},
     {"cal_ncolumns", &fNColumns, kInt},
+    {"cal_using_fadc", &fUsingFADC, kInt, 0, 1},
+    {"cal_ped_sample_low", &fPedSampLow, kInt, 0, 1},
+    {"cal_ped_sample_high", &fPedSampHigh, kInt, 0, 1},
+    {"cal_data_sample_low", &fDataSampLow, kInt, 0, 1},
+    {"cal_data_sample_high", &fDataSampHigh, kInt, 0, 1},
     {0}
   };
   gHcParms->LoadParmValues((DBRequest*)&list, prefix);
@@ -89,7 +100,9 @@ Int_t THcShowerArray::ReadDatabase( const TDatime& date )
   // Here read the 2-D arras of pedestals, gains, etc.
 
 
+  // Event by event amplitude and pedestal
   fA = new Double_t[fNelem];
+  fP = new Double_t[fNelem];
 
   return kOK;
 }
@@ -106,6 +119,7 @@ Int_t THcShowerArray::DefineVariables( EMode mode )
   RVarDef vars[] = {
     {"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"},
     {"a", "Raw ADC Amplitude", "fA"},
+    {"p", "Dynamic ADC Pedestal", "fP"},
     { 0 }
   };
 
@@ -173,7 +187,13 @@ Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
     THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
 
     // Should probably check that counter # is in range
-    fA[hit->fCounter-1] = hit->GetData(0);
+    if(fUsingFADC) {
+      fA[hit->fCounter-1] = hit->GetData(0,fPedSampLow,fPedSampHigh,
+					 fDataSampLow,fDataSampHigh);
+      fP[hit->fCounter-1] = hit->GetPedestal(0,fPedSampLow,fPedSampHigh);
+    } else {
+          fA[hit->fCounter-1] = hit->GetData(0);
+    }
 
     // Do other stuff like comparison to thresholds, signal hits, energy sums
 
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
index be52dad..681e886 100644
--- a/src/THcShowerArray.h
+++ b/src/THcShowerArray.h
@@ -55,6 +55,7 @@ public:
 protected:
 
   Double_t* fA;                 // [fNelem] ADC amplitude of blocks
+  Double_t* fP;                 // [fNelem] Event by event pedestals
 
   TClonesArray* fADCHits;	// List of ADC hits
 
@@ -64,6 +65,11 @@ protected:
   // Parameters
   Int_t fNRows;
   Int_t fNColumns;
+  Int_t fUsingFADC;		// != 0 if using FADC in sample mode
+  Int_t fPedSampLow;		// Sample range for
+  Int_t fPedSampHigh;		// dynamic pedestal
+  Int_t fDataSampLow;		// Sample range for
+  Int_t fDataSampHigh;		// sample integration
 
   Int_t fLayerNum;		// 2 for SHMS
   // Accumulators for pedestals go here
-- 
GitLab