From 1e6b00209e72fff3fe15584465f2f0484e7899d9 Mon Sep 17 00:00:00 2001
From: "Stephen A. Wood" <saw@jlab.org>
Date: Mon, 12 Oct 2015 10:38:08 -0400
Subject: [PATCH] Very basic support for the SHMS shower counter   Added as an
 optional sub detector THcShowerArray to the THcShower class   Get hit data
 with GetData method instead of hit member variables     so as to be
 compatible with FADC compatible hit list

---
 Makefile               |   1 +
 SConscript.py          |   1 +
 src/HallC_LinkDef.h    |   1 +
 src/THcRawShowerHit.h  |   1 +
 src/THcShower.cxx      |  49 +++++++--
 src/THcShower.h        |   5 +
 src/THcShowerArray.cxx | 218 +++++++++++++++++++++++++++++++++++++++++
 src/THcShowerArray.h   |  77 +++++++++++++++
 8 files changed, 347 insertions(+), 6 deletions(-)
 create mode 100644 src/THcShowerArray.cxx
 create mode 100644 src/THcShowerArray.h

diff --git a/Makefile b/Makefile
index 83f6afe..e69a57f 100644
--- a/Makefile
+++ b/Makefile
@@ -21,6 +21,7 @@ SRC  =  src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
 	src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \
 	src/THcSpacePoint.cxx src/THcDCTrack.cxx \
 	src/THcShower.cxx src/THcShowerPlane.cxx \
+	src/THcShowerArray.cxx \
 	src/THcRawShowerHit.cxx \
 	src/THcAerogel.cxx src/THcAerogelHit.cxx \
 	src/THcCherenkov.cxx src/THcCherenkovHit.cxx \
diff --git a/SConscript.py b/SConscript.py
index dc15d2a..7913179 100644
--- a/SConscript.py
+++ b/SConscript.py
@@ -17,6 +17,7 @@ hcheaders = Split("""
         src/THcDC.h src/THcDriftChamberPlane.h 
 	src/THcDriftChamber.h src/THcRawDCHit.h src/THcDCHit.h src/THcDCWire.h src/THcSpacePoint.h 
 	src/THcDCLookupTTDConv.h src/THcDCTimeToDistConv.h src/THcShower.h src/THcShowerPlane.h 
+        src/THcShowerArray.h
 	src/THcRawShowerHit.h src/THcAerogel.h src/THcAerogelHit.h src/THcCherenkov.h src/THcCherenkovHit.h
         src/THcGlobals.h src/THcDCTrack.h src/THcFormula.h
         src/THcRaster.h src/THcRasteredBeam.h src/THcRasterRawHit.h src/THcScalerEvtHandler.h
diff --git a/src/HallC_LinkDef.h b/src/HallC_LinkDef.h
index 93be610..d1efa4f 100644
--- a/src/HallC_LinkDef.h
+++ b/src/HallC_LinkDef.h
@@ -43,6 +43,7 @@
 #pragma link C++ class THcDCTrack+;
 #pragma link C++ class THcShower+;
 #pragma link C++ class THcShowerPlane+;
+#pragma link C++ class THcShowerArray+;
 #pragma link C++ class THcRawShowerHit+;
 #pragma link C++ class THcAerogel+;
 #pragma link C++ class THcAerogelHit+;
diff --git a/src/THcRawShowerHit.h b/src/THcRawShowerHit.h
index 76f0d14..934dced 100644
--- a/src/THcRawShowerHit.h
+++ b/src/THcRawShowerHit.h
@@ -7,6 +7,7 @@ class THcRawShowerHit : public THcRawHit {
 
  public:
   friend class THcShowerPlane;
+  friend class THcShowerArray;
 
   THcRawShowerHit(Int_t plane=0, Int_t counter=0) : THcRawHit(plane, counter), 
     fADC_pos(-1), fADC_neg(-1){
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index b547562..aae4120 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -39,6 +39,8 @@ THcShower::THcShower( const char* name, const char* description,
 {
   // Constructor
   fNLayers = 0;			// No layers until we make them
+  fNTotLayers = 0;
+  fHasArray = 0;
 
   fClusterList = new THcShowerClusterList;
 }
@@ -58,26 +60,29 @@ void THcShower::Setup(const char* name, const char* description)
   prefix[0] = tolower(GetApparatus()->GetName()[0]);
   prefix[1] = '\0';
 
+  
   string layernamelist;
+  fHasArray = 0;		// Flag for presence of fly's eye array
   DBRequest list[]={
     {"cal_num_layers", &fNLayers, kInt},
     {"cal_layer_names", &layernamelist, kString},
+    {"cal_array",&fHasArray, kInt,0, 1}, 
     {0}
   };
 
   gHcParms->LoadParmValues((DBRequest*)&list,prefix);
-
+  fNTotLayers = (fNLayers+(fHasArray!=0?1:0));
   vector<string> layer_names = vsplit(layernamelist);
 
-  if(layer_names.size() != (UInt_t) fNLayers) {
-    cout << "THcShower::Setup ERROR: Number of layers " << fNLayers 
+  if(layer_names.size() != fNTotLayers) {
+    cout << "THcShower::Setup ERROR: Number of layers " << fNTotLayers 
 	 << " doesn't agree with number of layer names "
 	 << layer_names.size() << endl;
     // Should quit.  Is there an official way to quit?
   }
 
-  fLayerNames = new char* [fNLayers];
-  for(UInt_t i=0;i<fNLayers;i++) {
+  fLayerNames = new char* [fNTotLayers];
+  for(UInt_t i=0;i<fNTotLayers;i++) {
     fLayerNames[i] = new char[layer_names[i].length()+1];
     strcpy(fLayerNames[i], layer_names[i].c_str());
   }
@@ -93,6 +98,16 @@ void THcShower::Setup(const char* name, const char* description)
     fPlanes[i] = new THcShowerPlane(fLayerNames[i], desc, i+1, this); 
 
   }
+  if(fHasArray) {
+    strcpy(desc, description);
+    strcat(desc, " Array ");
+    strcat(desc, fLayerNames[fNTotLayers-1]);
+
+    fArray = new THcShowerArray(fLayerNames[fNTotLayers-1], desc, fNTotLayers, this);
+  } else {
+    fArray = 0;
+  }
+
   delete [] desc;
 
   cout << "---------------------------------------------------------------\n";
@@ -127,6 +142,11 @@ THaAnalysisObject::EStatus THcShower::Init( const TDatime& date )
       return fStatus=status;
     }
   }
+  if(fHasArray) {
+    if((status = fArray->Init( date ))) {
+      return fStatus = status;
+    }
+  }
 
   char EngineDID[] = " CAL";
   EngineDID[0] = toupper(GetApparatus()->GetName()[0]);
@@ -171,9 +191,11 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
   prefix[0]=tolower(GetApparatus()->GetName()[0]);
   prefix[1]='\0';
 
+  fNegCols = fNLayers;		// If not defined assume tube on each end
+                                // for every layer
   {
     DBRequest list[]={
-      {"cal_num_neg_columns", &fNegCols, kInt},
+      {"cal_num_neg_columns", &fNegCols, kInt, 0, 1},
       {"cal_slop", &fSlop, kDouble},
       {"cal_fv_test", &fvTest, kInt,0,1},
       {"cal_fv_delta", &fvDelta, kDouble},
@@ -518,6 +540,9 @@ void THcShower::Clear(Option_t* opt)
   for(UInt_t ip=0;ip<fNLayers;ip++) {
     fPlanes[ip]->Clear(opt);
   }
+  if(fHasArray) {
+    fArray->Clear(opt);
+  }
 
   fNhits = 0;
   fNclust = 0;
@@ -552,6 +577,9 @@ Int_t THcShower::Decode( const THaEvData& evdata )
     for(UInt_t ip=0;ip<fNLayers;ip++) {
       nexthit = fPlanes[ip]->AccumulatePedestals(fRawHitList, nexthit);
     }
+    if(fHasArray) {
+      nexthit = fArray->AccumulatePedestals(fRawHitList, nexthit);
+    }
     fAnalyzePedestals = 1;	// Analyze pedestals first normal events
     return(0);
   }
@@ -560,6 +588,9 @@ Int_t THcShower::Decode( const THaEvData& evdata )
     for(UInt_t ip=0;ip<fNLayers;ip++) {
       fPlanes[ip]->CalculatePedestals();
     }
+    if(fHasArray) {
+      fArray->CalculatePedestals();
+    }
     fAnalyzePedestals = 0;	// Don't analyze pedestals next event
   }
 
@@ -568,6 +599,10 @@ Int_t THcShower::Decode( const THaEvData& evdata )
     nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
     fEtot += fPlanes[ip]->GetEplane();
   }
+  if(fHasArray) {
+    nexthit = fArray->ProcessHits(fRawHitList, nexthit);
+    fEtot += fArray->GetEplane();
+  }
   THcHallCSpectrometer *app = static_cast<THcHallCSpectrometer*>(GetApparatus());
   fEtotNorm=fEtot/(app->GetPcentral());
  
@@ -588,6 +623,8 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
 
   // Fill set of unclustered hits.
 
+  // Ignore shower array (SHMS) for now
+
   THcShowerHitSet HitSet;
 
   for(UInt_t j=0; j < fNLayers; j++) {
diff --git a/src/THcShower.h b/src/THcShower.h
index b2dddb9..fd5bac7 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -11,6 +11,7 @@
 #include "THaNonTrackingDetector.h"
 #include "THcHitList.h"
 #include "THcShowerPlane.h"
+#include "THcShowerArray.h"
 #include "TMath.h"
 
 
@@ -256,7 +257,10 @@ protected:
 
   char** fLayerNames;
   UInt_t fNLayers;	        // Number of layers in the calorimeter
+  UInt_t fNTotLayers;	        // Number of layers including array
+  UInt_t fHasArray;		// If !=0 fly's eye array behind preshower
   Double_t* fNLayerZPos;	// Z positions of fronts of layers
+  // Following apply to just sideways readout layers
   Double_t* BlockThick;		// Thickness of blocks
   UInt_t* fNBlocks;              // [fNLayers] number of blocks per layer
   UInt_t fNtotBlocks;            // Total number of shower counter blocks
@@ -287,6 +291,7 @@ protected:
   Double_t fDcor[2];
 
   THcShowerPlane** fPlanes;     // [fNLayers] Shower Plane objects
+  THcShowerArray* fArray;
 
   TClonesArray*  fTrackProj;    // projection of track onto plane
 
diff --git a/src/THcShowerArray.cxx b/src/THcShowerArray.cxx
new file mode 100644
index 0000000..7b684e6
--- /dev/null
+++ b/src/THcShowerArray.cxx
@@ -0,0 +1,218 @@
+//*-- Author :
+
+//////////////////////////////////////////////////////////////////////////
+//
+// THcShowerArray
+//
+//////////////////////////////////////////////////////////////////////////
+
+#include "THcShowerArray.h"
+#include "TClonesArray.h"
+#include "THcSignalHit.h"
+#include "THcGlobals.h"
+#include "THcParmList.h"
+#include "THcHitList.h"
+#include "THcShower.h"
+#include "THcRawShowerHit.h"
+#include "TClass.h"
+#include "math.h"
+#include "THaTrack.h"
+#include "THaTrackProj.h"
+
+#include <cstring>
+#include <cstdio>
+#include <cstdlib>
+#include <iostream>
+
+#include <fstream>
+using namespace std;
+
+ClassImp(THcShowerArray)
+
+//______________________________________________________________________________
+THcShowerArray::THcShowerArray( const char* name, 
+				const char* description,
+				const Int_t layernum,
+				THaDetectorBase* parent )
+  : THaSubDetector(name,description,parent)
+{
+  fADCHits = new TClonesArray("THcSignalHit",100);
+  fLayerNum = layernum;
+}
+
+THcShowerArray::~THcShowerArray()
+{
+  // Destructor
+  delete fADCHits;
+
+  delete [] fA;
+}
+
+//_____________________________________________________________________________
+THaAnalysisObject::EStatus THcShowerArray::Init( const TDatime& date )
+{
+  // Extra initialization for shower layer: set up DataDest map
+
+  if( IsZombie())
+    return fStatus = kInitError;
+
+  // How to get information for parent
+  //  if( GetParent() )
+  //    fOrigin = GetParent()->GetOrigin();
+
+  EStatus status;
+  if( (status=THaSubDetector::Init( date )) )
+    return fStatus = status;
+
+  return fStatus = kOK;
+
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::ReadDatabase( const TDatime& date )
+{
+
+  char prefix[2];
+  prefix[0]=tolower(GetParent()->GetPrefix()[0]);
+  prefix[1]='\0';
+
+  cout << "Parent name: " << GetParent()->GetPrefix() << endl;
+  DBRequest list[]={
+    {"cal_nrows", &fNRows, kInt},
+    {"cal_ncolumns", &fNColumns, kInt},
+    {0}
+  };
+  gHcParms->LoadParmValues((DBRequest*)&list, prefix);
+
+  fNelem = fNRows*fNColumns;
+
+  // Here read the 2-D arras of pedestals, gains, etc.
+
+
+  fA = new Double_t[fNelem];
+
+  return kOK;
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::DefineVariables( EMode mode )
+{
+  // Initialize global variables
+
+  if( mode == kDefine && fIsSetup ) return kOK;
+  fIsSetup = ( mode == kDefine );
+
+  // Register variables in global list
+  RVarDef vars[] = {
+    {"adchits", "List of ADC hits", "fADCHits.THcSignalHit.GetPaddleNumber()"},
+    {"a", "Raw ADC Amplitude", "fA"},
+    { 0 }
+  };
+
+  return DefineVarsFromList( vars, mode );
+}
+
+//_____________________________________________________________________________
+void THcShowerArray::Clear( Option_t* )
+{
+  // Clears the hit lists
+  fADCHits->Clear();
+
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::Decode( const THaEvData& evdata )
+{
+  // Doesn't actually get called.  Use Fill method instead
+
+  return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::CoarseProcess( TClonesArray& tracks )
+{
+
+  // Nothing is done here. See ProcessHits method instead.
+  //  
+
+ return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::FineProcess( TClonesArray& tracks )
+{
+  return 0;
+}
+
+//_____________________________________________________________________________
+Int_t THcShowerArray::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
+{
+  // Extract the data for this layer from hit list
+  // Assumes that the hit list is sorted by layer, so we stop when the
+  // plane doesn't agree and return the index for the next hit.
+
+  THcShower* fParent;
+  fParent = (THcShower*) GetParent();
+
+  // Initialize variables.
+
+  fADCHits->Clear();
+
+  for(Int_t i=0;i<fNelem;i++) {
+    fA[i] = 0;
+  }
+
+  // Process raw hits. Get ADC hits for the plane, assign variables for each
+  // channel.
+
+  Int_t nrawhits = rawhits->GetLast()+1;
+
+  Int_t ihit = nexthit;
+
+  while(ihit < nrawhits) {
+    THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
+
+    // Should probably check that counter # is in range
+    fA[hit->fCounter-1] = hit->fADC_pos;
+
+    // Do other stuff like comparison to thresholds, signal hits, energy sums
+
+    ihit++;
+  }
+
+  return(ihit);
+}
+//_____________________________________________________________________________
+Int_t THcShowerArray::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
+{
+  // Doesn't do anything yet except skip over hits
+
+  Int_t nrawhits = rawhits->GetLast()+1;
+
+  Int_t ihit = nexthit;
+
+  while(ihit < nrawhits) {
+    THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
+
+    // OK for hit list sorted by layer.
+    if(hit->fPlane > fLayerNum) {
+      break;
+    }
+    ihit++;
+  }
+  fNPedestalEvents++;
+
+  return(ihit);
+}
+//_____________________________________________________________________________
+void THcShowerArray::CalculatePedestals( )
+{
+  // Doesn't do anything yet
+
+}
+//_____________________________________________________________________________
+void THcShowerArray::InitializePedestals( )
+{
+  // Doesn't do anything yet
+  fNPedestalEvents = 0;
+} 
diff --git a/src/THcShowerArray.h b/src/THcShowerArray.h
new file mode 100644
index 0000000..be52dad
--- /dev/null
+++ b/src/THcShowerArray.h
@@ -0,0 +1,77 @@
+#ifndef ROOT_THcShowerArray
+#define ROOT_THcShowerArray 
+
+//////////////////////////////////////////////////////////////////////////////
+//                         
+// THcShowerArray
+//
+// A Hall C Fly's Eye Shower Array
+//
+// Subdetector for the fly's eye part of the SHMS shower counter.
+// 
+//////////////////////////////////////////////////////////////////////////////
+
+#include "THaSubDetector.h"
+#include "TClonesArray.h"
+
+#include <iostream>
+
+#include <fstream>
+
+class THaEvData;
+class THaSignalHit;
+
+class THcShowerArray : public THaSubDetector {
+
+public:
+  THcShowerArray( const char* name, const char* description,
+		  Int_t planenum, THaDetectorBase* parent = NULL);
+  virtual ~THcShowerArray();
+
+  virtual void    Clear( Option_t* opt="" );
+  virtual Int_t Decode( const THaEvData& );
+  virtual EStatus Init( const TDatime& run_time );
+
+  virtual Int_t CoarseProcess( TClonesArray& tracks );
+  virtual Int_t FineProcess( TClonesArray& tracks );
+  Bool_t   IsTracking() { return kFALSE; }
+  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( );
+  virtual void  InitializePedestals( );
+
+  //  Double_t fSpacing;   not used
+
+  TClonesArray* fParentHitList;
+
+  // Return zero for now
+  Double_t GetEplane() {
+    return 0.0;
+  };
+
+
+protected:
+
+  Double_t* fA;                 // [fNelem] ADC amplitude of blocks
+
+  TClonesArray* fADCHits;	// List of ADC hits
+
+  Int_t fNPedestalEvents;	/* Pedestal event counter */
+  Int_t fMinPeds;		/* Only analyze/update if num events > */
+
+  // Parameters
+  Int_t fNRows;
+  Int_t fNColumns;
+
+  Int_t fLayerNum;		// 2 for SHMS
+  // Accumulators for pedestals go here
+  // 2D arrays
+
+  virtual Int_t  ReadDatabase( const TDatime& date );
+  virtual Int_t  DefineVariables( EMode mode = kDefine );
+  //virtual void  InitializePedestals( );
+  ClassDef(THcShowerArray,0); // Fly;s Eye calorimeter array
+};
+#endif
-- 
GitLab