Pedestal work.

Use cut file to determine what is a pedestal event.  Pedestal events
only proceed to Decode.
In Decode (Just THcHodoscope for now) loop through each plane and accumulate
pedestal data.  Calculate the pedestals at the first non-pedestal event.
Pedestals are just calculated, but nothing is done with them.

examples/hodtest.C

@@ -64,7 +64,7 @@
   analyzer->SetEvent( event );
   analyzer->SetOutFile( "hodtest.root" );
   analyzer->SetOdefFile("output.def");
-  //  analyzer->SetCutFile("cuts_example.def");        // optional
+  analyzer->SetCutFile("hodtest_cuts.def");        // optional
   
   // File to record cuts accounting information
   //  analyzer->SetSummaryFile("summary_example.log"); // optional

examples/hodtest_cuts.def

+# Demo cuts for hodtest
+#
+
+Block: RawDecode
+
+Pedestal_event    g.evtyp==4
+RawDecode_master  1
+
+Block: Decode
+Decode_master     !xPedestal_event
+

src/THcAnalyzer.h

@@ -16,7 +16,11 @@ public:
   THcAnalyzer();
   virtual ~THcAnalyzer();
 
+  void SetPedestalEvtype( Int_t evtype ) { fPedestalEvtype = evtype; }
+
 protected:
+
+  Int_t fPedestalEvtype;
     
 private:
   //  THcAnalyzer( const THcAnalyzer& );

src/THcHodoscope.cxx

@@ -13,6 +13,8 @@
 #include "THaEvData.h"
 #include "THaDetMap.h"
 #include "THcDetectorMap.h"
+#include "THaGlobals.h"
+#include "THaCutList.h"
 #include "THcGlobals.h"
 #include "THcParmList.h"
 #include "VarDef.h"
@@ -330,6 +332,21 @@ Int_t THcHodoscope::Decode( const THaEvData& evdata )
   // Get the Hall C style hitlist (fRawHitList) for this event
   Int_t nhits = THcHitList::DecodeToHitList(evdata);
 
+  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;
   for(Int_t ip=0;ip<fNPlanes;ip++) {

src/THcHodoscope.h

@@ -41,6 +41,8 @@ public:
   THcHodoscope();  // for ROOT I/O
 protected:
 
+  Int_t fAnalyzePedestals;
+
   // Calibration
 
   // Per-event data

src/THcScintillatorPlane.cxx

@@ -87,7 +87,7 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
 
   // See what file it looks for
   
-  static const char* const here = "ReadDatabase()";
+  //  static const char* const here = "ReadDatabase()";
   char prefix[2];
   char parname[100];
   
@@ -124,7 +124,7 @@ Int_t THcScintillatorPlane::ReadDatabase( const TDatime& date )
 
 
   // Create arrays to hold results here
-
+  InitializePedestals();
 
   return kOK;
 }
@@ -187,6 +187,8 @@ Int_t THcScintillatorPlane::FineProcess( TClonesArray& tracks )
 {
   return 0;
 }
+
+//_____________________________________________________________________________
 Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
 {
   // Extract the data for this plane from hit list
@@ -271,6 +273,98 @@ Int_t THcScintillatorPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
   return(ihit);
 }
 
+//_____________________________________________________________________________
+Int_t THcScintillatorPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
+{
+  // Extract the data for this plane from hit list, accumulating into
+  // arrays for calculating pedestals.
+
+  Int_t nrawhits = rawhits->GetLast()+1;
+  // cout << "THcScintillatorPlane::AcculatePedestals " << fPlaneNum << " " << nexthit << "/" << nrawhits << endl;
+
+  Int_t ihit = nexthit;
+  while(ihit < nrawhits) {
+    THcHodoscopeHit* hit = (THcHodoscopeHit *) rawhits->At(ihit);
+    if(hit->fPlane > fPlaneNum) {
+      break;
+    }
+    Int_t element = hit->fCounter - 1; // Should check if in range
+    Int_t adcpos = hit->fADC_pos;
+    Int_t adcneg = hit->fADC_neg;
+
+    if(adcpos <= fPosPedLimit[element]) {
+      fPosPedSum[element] += adcpos;
+      fPosPedSum2[element] += adcpos*adcpos;
+      fPosPedCount[element]++;
+      if(fPosPedCount[element] == fMinPeds/5) {
+	fPosPedLimit[element] = 100 + fPosPedSum[element]/fPosPedCount[element];
+      }
+    }
+    if(adcneg <= fNegPedLimit[element]) {
+      fNegPedSum[element] += adcneg;
+      fNegPedSum2[element] += adcneg*adcneg;
+      fNegPedCount[element]++;
+      if(fNegPedCount[element] == fMinPeds/5) {
+	fNegPedLimit[element] = 100 + fNegPedSum[element]/fNegPedCount[element];
+      }
+    }
+    ihit++;
+  }
+
+  fNPedestalEvents++;
+
+  return(ihit);
+}
+
+//_____________________________________________________________________________
+void THcScintillatorPlane::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++) {
     
+    // Positive tubes
+    fPosPed[i] = ((Double_t) fPosPedSum[i]) / TMath::Max(1, fPosPedCount[i]);
+    fPosThresh[i] = fPosPed[i] + 15;
+
+    // Negative tubes
+    fNegPed[i] = ((Double_t) fNegPedSum[i]) / TMath::Max(1, fNegPedCount[i]);
+    fNegThresh[i] = fNegPed[i] + 15;
+
+    //    cout << i+1 << " " << fPosPed[i] << " " << fNegPed[i] << endl;
+  }
+  //  cout << " " << endl;
   
-  
+}
+
+//_____________________________________________________________________________
+void THcScintillatorPlane::InitializePedestals( )
+{
+  fNPedestalEvents = 0;
+  fMinPeds = 500; 		// In engine, this is set in parameter file
+  fPosPedSum = new Int_t [fNelem];
+  fPosPedSum2 = new Int_t [fNelem];
+  fPosPedLimit = new Int_t [fNelem];
+  fPosPedCount = new Int_t [fNelem];
+  fNegPedSum = new Int_t [fNelem];
+  fNegPedSum2 = new Int_t [fNelem];
+  fNegPedLimit = new Int_t [fNelem];
+  fNegPedCount = new Int_t [fNelem];
+
+  fPosPed = new Double_t [fNelem];
+  fNegPed = new Double_t [fNelem];
+  fPosThresh = new Double_t [fNelem];
+  fNegThresh = new Double_t [fNelem];
+  for(Int_t i=0;i<fNelem;i++) {
+    fPosPedSum[i] = 0;
+    fPosPedSum2[i] = 0;
+    fPosPedLimit[i] = 1000;	// In engine, this are set in parameter file
+    fPosPedCount[i] = 0;
+    fNegPedSum[i] = 0;
+    fNegPedSum2[i] = 0;
+    fNegPedLimit[i] = 1000;	// In engine, this are set in parameter file
+    fNegPedCount[i] = 0;
+  }
+}
+

src/THcScintillatorPlane.h

@@ -35,6 +35,8 @@ class THcScintillatorPlane : public THaSubDetector {
   virtual Bool_t   IsPid()      { return kFALSE; }
 
   virtual Int_t ProcessHits(TClonesArray* rawhits, Int_t nexthit);
+  virtual Int_t AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit);
+  virtual void  CalculatePedestals( );
 
   Double_t fSpacing;
 
@@ -51,11 +53,31 @@ class THcScintillatorPlane : public THaSubDetector {
   TClass* fNegTDCHitsClass;
   TClass* fPosADCHitsClass;
   TClass* fNegADCHitsClass;
-
-  Int_t fPlaneNum;
-
+  
+  Int_t fPlaneNum;		/* Which plane am I 1-4 */
+  Int_t fNelem;			/* Need since we don't inherit from 
+				 detector base class */
+  Int_t fNPedestalEvents;	/* Number of pedestal events */
+  Int_t fMinPeds;		/* Only analyze/update if num events > */
+  Int_t *fPosPedSum;		/* Accumulators for pedestals */
+  Int_t *fPosPedSum2;
+  Int_t *fPosPedLimit;
+  Int_t *fPosPedCount;
+  Int_t *fNegPedSum;
+  Int_t *fNegPedSum2;
+  Int_t *fNegPedLimit;
+  Int_t *fNegPedCount;
+
+  Double_t *fPosPed;
+  Double_t *fPosSig;
+  Double_t *fPosThresh;
+  Double_t *fNegPed;
+  Double_t *fNegSig;
+  Double_t *fNegThresh;
+  
   virtual Int_t  ReadDatabase( const TDatime& date );
   virtual Int_t  DefineVariables( EMode mode = kDefine );
+  virtual void  InitializePedestals( );
 
   ClassDef(THcScintillatorPlane,0)
 };