Added clustering of the HMS calorimeter hits

Makefile

@@ -21,7 +21,7 @@ SRC  =  src/THcInterface.cxx src/THcParmList.cxx src/THcAnalyzer.cxx \
 	src/THcDCWire.cxx \
 	src/THcDCLookupTTDConv.cxx src/THcDCTimeToDistConv.cxx \
 	src/THcShower.cxx src/THcShowerPlane.cxx \
-	src/THcShowerHit.cxx \
+	src/THcRawShowerHit.cxx \
 	src/THcAerogel.cxx src/THcAerogelHit.cxx
 
 # Name of your package. 

src/HallC_LinkDef.h

@@ -28,7 +28,7 @@
 #pragma link C++ class THcDCTimeToDistConv+;
 #pragma link C++ class THcShower+;
 #pragma link C++ class THcShowerPlane+;
-#pragma link C++ class THcShowerHit+;
+#pragma link C++ class THcRawShowerHit+;
 #pragma link C++ class THcAerogel+;
 #pragma link C++ class THcAerogelHit+;
 

src/THcShower.cxx

@@ -10,6 +10,7 @@
 ///////////////////////////////////////////////////////////////////////////////
 
 #include "THcShower.h"
+#include "THcShowerCluster.h"
 #include "THaEvData.h"
 #include "THaDetMap.h"
 #include "THcDetectorMap.h"
@@ -105,7 +106,7 @@ 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
 
-  THcHitList::InitHitList(fDetMap, "THcShowerHit", 100);
+  THcHitList::InitHitList(fDetMap, "THcRawShowerHit", 100);
 
   EStatus status;
   if( (status = THaNonTrackingDetector::Init( date )) )
@@ -186,6 +187,8 @@ Int_t THcShower::ReadDatabase( const TDatime& date )
     gHcParms->LoadParmValues((DBRequest*)&list, prefix);
   }
 
+  //Caution! Z positions (fronts) are off in hcal.param! Correct later on.
+
   YPos = new Double_t* [fNLayers];
   for(Int_t i=0;i<fNLayers;i++) {
     YPos[i] = new Double_t [fNBlocks[i]];
@@ -470,10 +473,10 @@ Int_t THcShower::Decode( const THaEvData& evdata )
     nexthit = fPlanes[ip]->ProcessHits(fRawHitList, nexthit);
   }
 
-  //   fRawHitList is TClones array of THcShowerHit objects
+  //   fRawHitList is TClones array of THcRawShowerHit objects
   //  cout << "THcShower::Decode: Shower raw hit list:" << endl;
   //  for(Int_t ihit = 0; ihit < fNRawHits ; ihit++) {
-  //    THcShowerHit* hit = (THcShowerHit *) fRawHitList->At(ihit);
+  //    THcRawShowerHit* hit = (THcRawShowerHit *) fRawHitList->At(ihit);
   //    cout << ihit << " : " << hit->fPlane << ":" << hit->fCounter << " : "
   //	 << hit->fADC_pos << " " << hit->fADC_neg << " "  << endl;
   //  }
@@ -508,8 +511,70 @@ Int_t THcShower::CoarseProcess( TClonesArray&  ) //tracks
   
   cout << "THcShower::CoarseProcess called ---------------------------" <<endl;
 
-  ApplyCorrections();
+  //  ApplyCorrections();
+
+  //
+  // Clustering of hits.
+  //
+
+  THcShowerHitList HitList;                    //list of unclusterd hits
+
+  for(Int_t j=0; j < fNLayers; j++) {
+
+   //cout << "Plane " << j << "  Eplane = " << fPlanes[j]->GetEplane() << endl;
 
+    for (Int_t i=0; i<fNBlocks[j]; i++) {
+
+      Float_t Edep = fPlanes[j]->GetEmean(i);
+      if (Edep > 0.) {                                    //hit
+	Float_t y = YPos[j][i] + BlockThick[j]/2.;        //top + thick/2
+	Float_t z = fNLayerZPos[j] + BlockThick[j]/2.;    //front + thick/2
+	THcShowerHit* hit = new THcShowerHit(i,j,y,z,Edep);
+
+	HitList.push_back(hit);
+
+	//cout << "Hit: Edep = " << Edep << " Y = " << y << " Z = " << z <<
+	//	  " Block " << i << " Layer " << j << endl;
+      };
+
+    }
+  }
+
+  //Print out hits before clustering.
+  //
+  cout << "Total hits:     " << HitList.size() << endl;
+  for (unsigned int i=0; i!=HitList.size(); i++) {
+    cout << "unclustered hit " << i << ": ";
+    (*(HitList.begin()+i))->show();
+  }
+
+  THcShowerClusterList* ClusterList = new THcShowerClusterList;
+  ClusterList->ClusterHits(HitList);
+
+  //Print out the cluster list.
+  //
+  cout << "Cluster_list size: " << (*ClusterList).NbClusters() << endl;
+
+  for (unsigned int i=0; i!=(*ClusterList).NbClusters(); i++) {
+
+    THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+
+    cout << "Cluster #" << i 
+         <<":  E=" << (*cluster).clE() 
+         << "  Epr=" << (*cluster).clEpr()
+         << "  Y=" << (*cluster).clY()
+         << "  Z=" << (*cluster).clZ()
+         << "  size=" << (*cluster).clSize()
+         << endl;
+
+    for (unsigned int j=0; j!=(*cluster).clSize(); j++) {
+      THcShowerHit* hit = (*cluster).ClusteredHit(j);
+      cout << "  hit #" << j << ":  "; (*hit).show();
+    }
+
+  }
+   
+  cout << "THcShower::CoarseProcess return ---------------------------" <<endl;
   return 0;
 }
 

src/THcShower.h

@@ -10,7 +10,6 @@
 #include "TClonesArray.h"
 #include "THaNonTrackingDetector.h"
 #include "THcHitList.h"
-#include "THcShowerHit.h"
 #include "THcShowerPlane.h"
 
 class THaScCalib;

src/THcShowerHit.cxx

-///////////////////////////////////////////////////////////////////////////////
-//                                                                           //
-// THcShowerHit                                                           //
-//                                                                           //
-// Class representing a single raw hit for a hodoscope paddle                //
-//                                                                           //
-// Contains plane, counter and pos/neg adc and tdc values                    //
-//                                                                           //
-///////////////////////////////////////////////////////////////////////////////
-
-#include "THcShowerHit.h"
-
-using namespace std;
-
-
-void THcShowerHit::SetData(Int_t signal, Int_t data) {
-  if(signal==0) {
-    fADC_pos = data;
-  } else if (signal==1) {
-    fADC_neg = data;
-  } 
-}
-
-Int_t THcShowerHit::GetData(Int_t signal) {
-  if(signal==0) {
-    return(fADC_pos);
-  } else if (signal==1) {
-    return(fADC_neg);
-  } 
-
-  return(-1); // Actually should throw exception
-}
-
-#if 0
-Int_t THcShowerHit::Compare(const TObject* obj) const
-{
-  // Compare to sort by plane and counter
-
-  const THcShowerHit* hit = dynamic_cast<const THcShowerHit*>(obj);
-
-  if(!hit) return -1;
-  Int_t p1 = fPlane;
-  Int_t p2 = hit->fPlane;
-  if(p1 < p2) return -1;
-  else if(p1 > p2) return 1;
-  else {
-    Int_t c1 = fCounter;
-    Int_t c2 = hit->fCounter;
-    if(c1 < c2) return -1;
-    else if (c1 == c2) return 0;
-    else return 1;
-  }
-}
-#endif
-//_____________________________________________________________________________
-THcShowerHit& THcShowerHit::operator=( const THcShowerHit& rhs )
-{
-  // Assignment operator.
-
-  THcRawHit::operator=(rhs);
-  if ( this != &rhs ) {
-    fPlane = rhs.fPlane;
-    fCounter = rhs.fCounter;
-    fADC_pos = rhs.fADC_pos;
-    fADC_neg = rhs.fADC_neg;
-
-  }
-  return *this;
-}
-
-
-//////////////////////////////////////////////////////////////////////////
-ClassImp(THcShowerHit)
-
- 

src/THcShowerHit.h

 #ifndef ROOT_THcShowerHit
 #define ROOT_THcShowerHit
 
-#include "THcRawHit.h"
+// HMS calorimeter hit, version 2
 
-class THcShowerHit : public THcRawHit {
+#include <vector>
+#include <iterator>
+#include <iostream>
+
+using namespace std;
+
+class THcShowerHit {           //HMS calorimeter hit class
+
+ private:
+  unsigned int fCol, fRow; //hit colomn and row
+  float fY, fZ;            //hit Y (vert.) and Z (along spect.axis) coordinates
+  float fE;                 //hit energy deposition
 
  public:
 
- THcShowerHit(Int_t plane=0, Int_t counter=0) : THcRawHit(plane, counter), 
-    fADC_pos(-1), fADC_neg(-1){
+  THcShowerHit() {             //default constructor
+  fCol=fRow=0;
+  fY=fZ=0.;
+  fE=0.;
   }
-  THcShowerHit& operator=( const THcShowerHit& );
-  virtual ~THcShowerHit() {}
 
-  virtual void Clear( Option_t* opt="" )
-    { fADC_pos = -1; fADC_neg = -1; }
+  THcShowerHit(unsigned int hRow, unsigned int hCol, float hY, float hZ,
+	       float hE) {
+    fRow=hRow;
+    fCol=hCol;
+    fY=hY;
+    fZ=hZ;
+    fE=hE;
+  }
 
-  void SetData(Int_t signal, Int_t data);
-  Int_t GetData(Int_t signal);
+  ~THcShowerHit() {
+    //    cout << " hit destructed" << endl;
+  }
 
-  //  virtual Bool_t  IsSortable () const {return kTRUE; }
-  //  virtual Int_t   Compare(const TObject* obj) const;
+  unsigned int hitColumn() {
+    return fCol;
+  }
 
-  Int_t fADC_pos;
-  Int_t fADC_neg;
+  unsigned int hitRow() {
+    return fRow;
+  }
 
- protected:
+  float hitY() {
+    return fY;
+  }
 
- private:
+  float hitZ() {
+    return fZ;
+  }
 
-  ClassDef(THcShowerHit, 0);	// Shower hit class
-};  
+  float hitE() {
+    return fE;
+  }
+
+  bool isNeighbour(THcShowerHit* hit1) {      //Is hit1 neighbouring this hit?
+    int dRow = fRow-(*hit1).fRow;
+    int dCol = fCol-(*hit1).fCol;
+    return abs(dRow)<2 && abs(dCol)<2;
+  }
+
+  //Print out hit information
+  //
+  void show() {
+    cout << "row=" << fRow << "  column=" << fCol << 
+      "  y=" << fY << "  z=" << fZ << "  E=" << fE << endl;
+  }
+
+};
+
+
+typedef vector<THcShowerHit*> THcShowerHitList;      //alias for hit container
+typedef THcShowerHitList::iterator THcShowerHitIt;   //and for its iterator
+
+//Purge sequence container of pointers. Found in Echel, v.2, p.253.
+//
+template<class Seq> void purge(Seq& c) {
+  typename Seq::iterator i;
+  for(i = c.begin(); i != c.end(); i++) {
+    delete *i;
+    *i = 0;
+  }
+}
 
 #endif
- 

src/THcShowerPlane.cxx

@@ -13,6 +13,7 @@
 #include "THcParmList.h"
 #include "THcHitList.h"
 #include "THcShower.h"
+#include "THcRawShowerHit.h"
 #include "TClass.h"
 #include "math.h"
 
@@ -253,7 +254,7 @@ Int_t THcShowerPlane::ProcessHits(TClonesArray* rawhits, Int_t nexthit)
   //cout << "nrawhits =  " << nrawhits << endl;
   //cout << "nexthit =  " << nexthit << endl;
   while(ihit < nrawhits) {
-    THcShowerHit* hit = (THcShowerHit *) rawhits->At(ihit);
+    THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
 
     //cout << "fplane =  " << hit->fPlane << " Num = " << fLayerNum << endl;
     if(hit->fPlane > fLayerNum) {
@@ -308,7 +309,7 @@ Int_t THcShowerPlane::AccumulatePedestals(TClonesArray* rawhits, Int_t nexthit)
   //  cout << "THcScintillatorPlane::AcculatePedestals " << fLayerNum << " " << nexthit << "/" << nrawhits << endl;
   Int_t ihit = nexthit;
   while(ihit < nrawhits) {
-    THcShowerHit* hit = (THcShowerHit *) rawhits->At(ihit);
+    THcRawShowerHit* hit = (THcRawShowerHit *) rawhits->At(ihit);
     //cout << "fPlane =  " << hit->fPlane << " Limit = " << fPlaneNum << endl;
     if(hit->fPlane > fLayerNum) {
       break;

src/THcShowerPlane.h

@@ -46,6 +46,14 @@ class THcShowerPlane : public THaSubDetector {
 
   TClonesArray* fParentHitList;
 
+  Float_t GetEplane() {
+    return fEplane;
+  };
+
+  Float_t GetEmean(Int_t i) {
+    return fEmean[i];
+  };
+
  protected:
 
   Float_t*   fA_Pos;         // [fNelem] Array of ADC amplitudes of blocks