diff --git a/Makefile b/Makefile
index 422e44b8456bc4d1847309026d5fb6da4dcf2d4b..ed806bd80b2423e2ce18a92833fcff2b00f9e792 100644
--- a/Makefile
+++ b/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.
diff --git a/src/HallC_LinkDef.h b/src/HallC_LinkDef.h
index ce520fa067dda33b10aebdad9bebadeb2336f719..dc2ada8fcb2ebed87ccdb92e044a8542bb616972 100644
--- a/src/HallC_LinkDef.h
+++ b/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+;
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index 7525a5f98682793df66e3a6dcd893b568be66a8e..db9b4ab11e929a4aeba7d97e7ff30bd6e7846171 100644
--- a/src/THcShower.cxx
+++ b/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;
}
diff --git a/src/THcShower.h b/src/THcShower.h
index 08ba824e3b0410e1e8a45222fa89d8063ef9bd15..dbab5b596f715f803dd0bdd12509db2e94377d5a 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -10,7 +10,6 @@
#include "TClonesArray.h"
#include "THaNonTrackingDetector.h"
#include "THcHitList.h"
-#include "THcShowerHit.h"
#include "THcShowerPlane.h"
class THaScCalib;
diff --git a/src/THcShowerHit.cxx b/src/THcShowerHit.cxx
deleted file mode 100644
index c27357235534b2b8c67d13cdbcba0984c6fb8b09..0000000000000000000000000000000000000000
--- a/src/THcShowerHit.cxx
+++ /dev/null
@@ -1,75 +0,0 @@
-///////////////////////////////////////////////////////////////////////////////
-// //
-// 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)
-
-
diff --git a/src/THcShowerHit.h b/src/THcShowerHit.h
index 7d1a5f08ee650cf8eb3a96dd82ce4a9ff414171f..d2504d23868767147fa1494df9cea82529ed4cbf 100644
--- a/src/THcShowerHit.h
+++ b/src/THcShowerHit.h
@@ -1,36 +1,89 @@
#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
-
diff --git a/src/THcShowerPlane.cxx b/src/THcShowerPlane.cxx
index bfe626910e8103c9fc22846df5fb3dcd2c68fa17..a6d642f3b8aa2506c3baa8fd85422826049c3d2d 100644
--- a/src/THcShowerPlane.cxx
+++ b/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;
diff --git a/src/THcShowerPlane.h b/src/THcShowerPlane.h
index 2e04130c42b6ffdb6531ee73e82dc221ac6a2bb4..149cb2014e7e6466d3618ef63a009270c799be7f 100644
--- a/src/THcShowerPlane.h
+++ b/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