From f1b25902d866ec04625b34d6af68a6f1b379815d Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Mon, 19 Jan 2015 11:30:01 +0400
Subject: [PATCH] Improve THcShower class structure and cleanup Clean up
THcShower.h and THcShower.cxx. typedef vector<THcShowerHit*>
THcShowerHitList by typedef set<THcShowerHit*> THcShowerHitList. Remove
inheritance of class THcShowerCluster from THcShowerHitList. Define
THcShowerCluster as THcShowerHitList, aka set<THcShowerHit*>. Remove
inheritance of class THcShowerClusterList from THcShClusterList,
a.k.a. vector<THcShowerCluster*> container. Define THcShowerClusterList as
vector<THcShowerCluster*>. Rename THcShowerHitSet by THcShowerHitList.
Rename HitList (related to the HMS calorimeter) to HitSet. Add operator< in
THcShowerHit class, in order to have THcShowerHitSet objects to be
properly sorted.
---
examples/PARAM/52949/hcal.param.vt.52949 | 2 +-
src/THcShower.cxx | 212 +++++++++++++++---
src/THcShower.h | 272 +++--------------------
3 files changed, 219 insertions(+), 267 deletions(-)
diff --git a/examples/PARAM/52949/hcal.param.vt.52949 b/examples/PARAM/52949/hcal.param.vt.52949
index 8e5f719..3970fb3 100644
--- a/examples/PARAM/52949/hcal.param.vt.52949
+++ b/examples/PARAM/52949/hcal.param.vt.52949
@@ -4,7 +4,7 @@ hcal_slop = 7.5
;Turn on HMS cal. fiducial volume cut. 0="no cut"
;Default hcal_fv_test=0
-hcal_fv_test = 0
+hcal_fv_test = 1
hcal_pos_cal_const =0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001
0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001,0.001
diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index bb89477..3ed60a6 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -29,6 +29,7 @@
#include <cstdio>
#include <cstdlib>
#include <iostream>
+#include <numeric>
using namespace std;
@@ -521,7 +522,14 @@ void THcShower::Clear(Option_t* opt)
fEtot = 0.;
fEtotNorm = 0.;
- fClusterList->purge();
+ // Purge cluster list
+
+ for (THcShowerClusterListIt i=fClusterList->begin(); i!=fClusterList->end();
+ ++i) {
+ delete *i;
+ *i = 0;
+ }
+ fClusterList->clear();
}
@@ -572,12 +580,12 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
//
// Apply corrections and reconstruct the complete hits.
- // Clustering of hits.
- //
+ // Clustering of hits.
+ //
- // Fill list of unclustered hits.
+ // Fill set of unclustered hits.
- THcShowerHitList HitList;
+ THcShowerHitSet HitSet;
for(UInt_t j=0; j < fNLayers; j++) {
@@ -606,13 +614,13 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
THcShowerHit* hit = new THcShowerHit(i,j,x,z,Edep,Epos,Eneg);
- HitList.push_back(hit);
+ HitSet.insert(hit); //<set> version
}
}
}
- fNhits = HitList.size();
+ fNhits = HitSet.size();
//Debug output, print out hits before clustering.
@@ -620,17 +628,18 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
cout << "---------------------------------------------------------------\n";
cout << "Debug output from THcShower::CoarseProcess\n";
cout << " List of unclustered hits. Total hits: " << fNhits << endl;
+ THcShowerHitIt it = HitSet.begin(); //<set> version
for (Int_t i=0; i!=fNhits; i++) {
cout << " hit " << i << ": ";
- (*(HitList.begin()+i))->show();
+ (*(it++))->show();
}
}
// Create list of clusters and fill it.
- fClusterList->ClusterHits(HitList);
+ ClusterHits(HitSet);
- fNclust = (*fClusterList).NbClusters(); //number of clusters
+ fNclust = (*fClusterList).size(); //number of clusters
//Debug output, print out the cluster list.
@@ -638,20 +647,23 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
cout << " Clustered hits. Number of clusters: " << fNclust << endl;
- for (Int_t i=0; i!=fNclust; i++) {
+ UInt_t i = 0;
+ for (THcShowerClusterListIt ppcl = (*fClusterList).begin();
+ ppcl != (*fClusterList).end(); ppcl++) {
- THcShowerCluster* cluster = (*fClusterList).ListedCluster(i);
- cout << " Cluster #" << i
- <<": E=" << (*cluster).clE()
- << " Epr=" << (*cluster).clEpr()
- << " X=" << (*cluster).clX()
- << " Z=" << (*cluster).clZ()
- << " size=" << (*cluster).clSize()
+ cout << " Cluster #" << i++
+ <<": E=" << clE(*ppcl)
+ << " Epr=" << clEpr(*ppcl)
+ << " X=" << clX(*ppcl)
+ << " Z=" << clZ(*ppcl)
+ << " size=" << (**ppcl).size()
<< endl;
- for (UInt_t j=0; j!=(*cluster).clSize(); j++) {
- THcShowerHit* hit = (*cluster).ClusteredHit(j);
- cout << " hit #" << j << ": "; (*hit).show();
+ Int_t j=0;
+ for (THcShowerClusterIt pph=(**ppcl).begin(); pph!=(**ppcl).end();
+ pph++) {
+ cout << " hit " << j++ << ": ";
+ (**pph).show();
}
}
@@ -664,6 +676,149 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
//-----------------------------------------------------------------------------
+void THcShower::ClusterHits(THcShowerHitSet& HitSet) {
+
+ // Collect hits from the HitSet into the clusters. The resultant clusters
+ // of hits are saved in the ClusterList.
+
+ while (HitSet.size() != 0) {
+
+ THcShowerCluster* cluster = new THcShowerCluster;
+
+ THcShowerHitIt it = HitSet.end();
+ (*cluster).insert(*(--it)); //Move the last hit from the hit list
+ HitSet.erase(it); //into the 1st cluster
+
+ bool clustered = true;
+
+ while (clustered) { //Proceed while a hit is clustered
+
+ clustered = false;
+
+ for (THcShowerHitIt i=HitSet.begin(); i!=HitSet.end(); ++i) {
+
+ for (THcShowerClusterIt k=(*cluster).begin(); k!=(*cluster).end();
+ k++) {
+
+ if ((**i).isNeighbour(*k)) {
+
+ (*cluster).insert(*i); //If the hit #i is neighbouring a hit
+ HitSet.erase(i); //in the cluster, then move it
+ //into the cluster.
+ clustered = true;
+ }
+
+ if (clustered) break;
+ } //k
+
+ if (clustered) break;
+ } //i
+
+ } //while clustered
+
+ fClusterList->push_back(cluster); //Put the cluster in the cluster list
+
+ } //While hit_list not exhausted
+
+};
+
+//-----------------------------------------------------------------------------
+
+Double_t addE(Double_t x, THcShowerHit* h) {
+ return x + h->hitE();
+}
+
+Double_t addX(Double_t x, THcShowerHit* h) {
+ return x + h->hitE() * h->hitX();
+}
+
+Double_t addZ(Double_t x, THcShowerHit* h) {
+ return x + h->hitE() * h->hitZ();
+}
+
+Double_t addEpr(Double_t x, THcShowerHit* h) {
+ return h->hitColumn() == 0 ? x + h->hitE() : x;
+}
+
+Double_t addEpos(Double_t x, THcShowerHit* h) {
+ return x + h->hitEpos();
+}
+
+Double_t addEneg(Double_t x, THcShowerHit* h) {
+ return x + h->hitEneg();
+}
+
+// X coordinate of center of gravity of cluster, calculated as hit energy
+// weighted average. Put X out of the calorimeter (-75 cm), if there is no
+// energy deposition in the cluster.
+//
+Double_t clX(THcShowerCluster* cluster) {
+ Double_t Etot = accumulate((*cluster).begin(),(*cluster).end(),0.,addE);
+ return (Etot != 0. ?
+ accumulate((*cluster).begin(),(*cluster).end(),0.,addX)/Etot : -75.);
+}
+
+// Z coordinate of center of gravity of cluster, calculated as a hit energy
+// weighted average. Put Z out of the calorimeter (0 cm), if there is no energy
+// deposition in the cluster.
+//
+Double_t clZ(THcShowerCluster* cluster) {
+ Double_t Etot = accumulate((*cluster).begin(),(*cluster).end(),0.,addE);
+ return (Etot != 0. ?
+ accumulate((*cluster).begin(),(*cluster).end(),0.,addZ)/Etot : 0.);
+}
+
+//Energy depostion in a cluster
+//
+Double_t clE(THcShowerCluster* cluster) {
+ return accumulate((*cluster).begin(),(*cluster).end(),0.,addE);
+}
+
+//Energy deposition in the Preshower (1st plane) for a cluster
+//
+Double_t clEpr(THcShowerCluster* cluster) {
+ return accumulate((*cluster).begin(),(*cluster).end(),0.,addEpr);
+}
+
+
+//Cluster energy deposition in plane iplane=0,..,3:
+// side=0 -- from positive PMTs only;
+// side=1 -- from negative PMTs only;
+// side=2 -- from postive and negative PMTs.
+//
+
+Double_t clEplane(THcShowerCluster* cluster, Int_t iplane, Int_t side) {
+
+ if (side!=0&&side!=1&&side!=2) {
+ cout << "*** Wrong Side in clEplane:" << side << " ***" << endl;
+ return -1;
+ }
+
+ THcShowerHitSet pcluster;
+ for (THcShowerHitIt it=(*cluster).begin(); it!=(*cluster).end(); ++it) {
+ if ((*it)->hitColumn() == iplane) pcluster.insert(*it);
+ }
+
+ Double_t Eplane;
+ switch (side) {
+ case 0 :
+ Eplane = accumulate(pcluster.begin(),pcluster.end(),0.,addEpos);
+ break;
+ case 1 :
+ Eplane = accumulate(pcluster.begin(),pcluster.end(),0.,addEneg);
+ break;
+ case 2 :
+ Eplane = accumulate(pcluster.begin(),pcluster.end(),0.,addE);
+ break;
+ default :
+ Eplane = 0. ;
+ }
+
+ return Eplane;
+}
+
+//-----------------------------------------------------------------------------
+
Int_t THcShower::MatchCluster(THaTrack* Track,
Double_t& XTrFront, Double_t& YTrFront)
{
@@ -719,16 +874,14 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
if (inFidVol) {
- // for (Int_t i=0; i<fNclust; i++) {
-
// Since hits and clusters are in reverse order (with respect to Engine),
// search backwards to be consistent with Engine.
//
for (Int_t i=fNclust-1; i>-1; i--) {
- THcShowerCluster* cluster = (*fClusterList).ListedCluster(i);
+ THcShowerCluster* cluster = *(fClusterList->begin()+i);
- Double_t dx = TMath::Abs( (*cluster).clX() - XTrFront );
+ Double_t dx = TMath::Abs( clX(cluster) - XTrFront );
if (dx <= (0.5*BlockThick[0] + fSlop)) {
fNtracks++; // number of shower tracks (Consistent with engine)
@@ -761,7 +914,7 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
cout << " Fiducial volume test: inFidVol = " << inFidVol << endl;
cout << " Matched cluster #" << mclust << ", delatX= " << deltaX
- << endl;
+ << endl;
cout << "---------------------------------------------------------------\n";
}
@@ -789,7 +942,8 @@ Float_t THcShower::GetShEnergy(THaTrack* Track) {
if (mclust >= 0) { // if there is a matched cluster
// Get matched cluster.
- THcShowerCluster* cluster = (*fClusterList).ListedCluster(mclust);
+
+ THcShowerCluster* cluster = *(fClusterList->begin()+mclust);
// Correct track energy depositions for the impact coordinate.
@@ -809,8 +963,8 @@ Float_t THcShower::GetShEnergy(THaTrack* Track) {
corneg = 0.;
}
- Etrk += (*cluster).clEplane(ip,0) * corpos;
- Etrk += (*cluster).clEplane(ip,1) * corneg;
+ Etrk += clEplane(cluster,ip,0) * corpos;
+ Etrk += clEplane(cluster,ip,1) * corneg;
} //over planes
diff --git a/src/THcShower.h b/src/THcShower.h
index 0e1880e..e52b424 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -16,7 +16,7 @@
// HMS calorimeter hits, version 2
-#include <vector>
+#include <set>
#include <iterator>
#include <iostream>
#include <memory>
@@ -86,12 +86,12 @@ public:
}
// Decide if a hit is neighbouring the current hit.
- // Two hits are neighbours if share a side or a corner.
+ // Two hits are neighbours if share a side or a corner,
+ // or in the same row but separated by no more than a block.
//
bool isNeighbour(THcShowerHit* hit1) { //Is hit1 neighbouring this hit?
Int_t dRow = fRow-(*hit1).fRow;
Int_t dCol = fCol-(*hit1).fCol;
- // return TMath::Abs(dRow)<2 && TMath::Abs(dCol)<2;
return (TMath::Abs(dRow)<2 && TMath::Abs(dCol)<2) ||
(dRow==0 && TMath::Abs(dCol)<3);
}
@@ -104,252 +104,36 @@ public:
<< " E=" << fE << " Epos=" << fEpos << " Eneg=" << fEneg << endl;
}
-};
-
-
-// Container (collection) of hits and its iterator.
-//
-typedef vector<THcShowerHit*> THcShowerHitList;
-typedef THcShowerHitList::iterator THcShowerHitIt;
-
-
-
-//HMS calorimeter hit cluster
-//
-class THcShowerCluster : THcShowerHitList {
-
- public:
-
- THcShowerCluster() {
- // cout << "Dummy THcShowerCluster object created" << endl;
- }
-
- ~THcShowerCluster() {
- for (THcShowerHitIt i = THcShowerHitList::begin();
- i != THcShowerHitList::end(); ++i) {
- delete *i;
- *i = 0;
- }
- }
-
- // Add a hit to the cluster hit list
- //
- void grow(THcShowerHit* hit) {
- THcShowerHitList::push_back(hit);
- }
-
- //Pointer to the hit #i in the cluster hit list
- //
- THcShowerHit* ClusteredHit(UInt_t i) {
- return * (THcShowerHitList::begin()+i);
- }
-
- //Print out a hit in the cluster
- //
- void showHit(UInt_t num) {
- (*(THcShowerHitList::begin()+num))->show();
- }
-
- // X coordinate of center of gravity of cluster,
- // calculated as hit energy weighted average.
- // Put X out of the calorimeter (-75 cm), if there is no energy deposition
- // in the cluster.
+ // Define < operator in order to fill in hit sets in a sorted manner.
//
- Double_t clX() {
- Double_t x_sum=0.;
- Double_t Etot=0.;
- for (THcShowerHitIt it=THcShowerHitList::begin();
- it!=THcShowerHitList::end(); ++it) {
- x_sum += (*it)->hitX() * (*it)->hitE();
- Etot += (*it)->hitE();
- }
- // cout << "x_sum=" << x_sum << " Etot=" << Etot << endl;
- return (Etot != 0. ? x_sum/Etot : -75.);
- }
-
- // Z coordinate of center of gravity of cluster,
- // calculated as a hit energy weighted average.
- // Put Z out of the calorimeter (0 cm), if there is no energy deposition
- // in the cluster.
- //
- Double_t clZ() {
- Double_t z_sum=0.;
- Double_t Etot=0.;
- for (THcShowerHitIt it=THcShowerHitList::begin();
- it!=THcShowerHitList::end(); ++it) {
- z_sum += (*it)->hitZ() * (*it)->hitE();
- Etot += (*it)->hitE();
- }
- // cout << "z_sum=" << z_sum << " Etot=" << Etot << endl;
- return (Etot != 0. ? z_sum/Etot : 0.);
- }
-
- //Energy depostion in a cluster
- //
- Double_t clE() {
- // cout << "In ECl:" << endl;
- Double_t Etot=0.;
- for (THcShowerHitIt it=THcShowerHitList::begin();
- it!=THcShowerHitList::end(); ++it) {
- Etot += (*it)->hitE();
- }
- return Etot;
- }
-
- //Energy deposition in the Preshower (1st plane) for a cluster
- //
- Double_t clEpr() {
- Double_t Epr=0.;
- for (THcShowerHitIt it=THcShowerHitList::begin();
- it!=THcShowerHitList::end(); ++it) {
- if ((*it)->hitColumn() == 0) {
- Epr += (*it)->hitE();
- }
- }
- return Epr;
- }
-
- //Cluster energy deposition in plane iplane=0,..,3:
- // side=0 -- from positive PMTs only;
- // side=1 -- from negative PMTs only;
- // side=2 -- from postive and negative PMTs.
- //
-
- Double_t clEplane(Int_t iplane, Int_t side) {
-
- if (side!=0&&side!=1&&side!=2) {
- cout << "*** Wrong Side in clEplane:" << side << " ***" << endl;
- return -1;
- }
-
- Double_t Eplane=0.;
- for (THcShowerHitIt it=THcShowerHitList::begin();
- it!=THcShowerHitList::end(); ++it) {
-
- if ((*it)->hitColumn() == iplane)
-
- switch (side) {
- case 0 : Eplane += (*it)->hitEpos();
- break;
- case 1 : Eplane += (*it)->hitEneg();
- break;
- case 2 : Eplane += (*it)->hitE();
- break;
- default : ;
- }
-
- }
-
- return Eplane;
- }
-
-
- //Cluster size (number of hits in the cluster).
- //
- UInt_t clSize() {
- return THcShowerHitList::size();
+ bool operator<(THcShowerHit rhs) const {
+ if (fCol != rhs.fCol)
+ return fCol < rhs.fCol;
+ else
+ return fRow < rhs.fRow;
}
};
-//-----------------------------------------------------------------------------
-//Alias for container of clusters and for its iterator
-//
-typedef vector<THcShowerCluster*> THcShClusterList;
-typedef THcShClusterList::iterator THcShClusterIt;
+//____________________________________________________________________________
-//List of clusters
+// Container (collection) of hits and its iterator.
//
-class THcShowerClusterList : private THcShClusterList {
-
- public:
-
- THcShowerClusterList() {
- // cout << "Dummy THcShowerClusterList object created" << endl;
- }
-
- ~THcShowerClusterList() {
- purge();
- }
+typedef set<THcShowerHit*> THcShowerHitSet;
+typedef THcShowerHitSet::iterator THcShowerHitIt;
- // Purge cluster list
- //
- void purge() {
- for (THcShClusterIt i = THcShClusterList::begin();
- i != THcShClusterList::end(); ++i) {
- delete *i;
- *i = 0;
- }
- THcShClusterList::clear();
- }
-
- //Put a cluster in the cluster list
- //
- void grow(THcShowerCluster* cluster) {
- THcShClusterList::push_back(cluster);
- }
-
- //Pointer to the cluster #i in the cluster list
- //
- THcShowerCluster* ListedCluster(UInt_t i) {
- return *(THcShClusterList::begin()+i);
- }
-
- //Cluster list size.
- //
- UInt_t NbClusters() {
- return THcShClusterList::size();
- }
+typedef THcShowerHitSet THcShowerCluster;
+typedef THcShowerCluster::iterator THcShowerClusterIt;
- //____________________________________________________________________________
-
- void ClusterHits(THcShowerHitList HitList) {
-
- // Collect hits from the HitList into the clusters. The resultant clusters
- // of hits are saved in the ClusterList.
-
- while (HitList.size() != 0) {
-
- THcShowerCluster* cluster = new THcShowerCluster;
-
- (*cluster).grow(*(HitList.end()-1)); //Move the last hit from the hit list
- HitList.erase(HitList.end()-1); //into the 1st cluster
- bool clustered = true;
-
- while (clustered) { //Proceed while a hit is clustered
-
- clustered = false;
-
- for (THcShowerHitIt i=HitList.begin(); i!=HitList.end(); ++i) {
-
- for (UInt_t k=0; k!=(*cluster).clSize(); k++) {
-
- if ((**i).isNeighbour((*cluster).ClusteredHit(k))) {
-
- (*cluster).grow(*i); //If the hit #i is neighbouring a hit
- HitList.erase(i); //in the cluster, then move it
- //into the cluster.
- clustered = true;
- }
-
- if (clustered) break;
- } //k
-
- if (clustered) break;
- } //i
-
- } //while clustered
-
- grow(cluster); //Put the cluster in the cluster list
-
- } //While hit_list not exhausted
-
- }
-
-};
+//______________________________________________________________________________
+//Alias for container of clusters and for its iterator
+//
+typedef vector<THcShowerCluster*> THcShowerClusterList;
+typedef THcShowerClusterList::iterator THcShowerClusterListIt;
+//______________________________________________________________________________
class THcShower : public THaNonTrackingDetector, public THcHitList {
@@ -513,6 +297,8 @@ protected:
// Cluster to track association method.
Int_t MatchCluster(THaTrack*, Double_t&, Double_t&);
+ void ClusterHits(THcShowerHitSet& HitSet);
+
friend class THcShowerPlane; //to access debug flags.
ClassDef(THcShower,0) // Shower counter detector
@@ -520,4 +306,16 @@ protected:
///////////////////////////////////////////////////////////////////////////////
+// Auxiliary methods to be used with the hit and cluster containers.
+
+Double_t addE(Double_t x, THcShowerHit* h);
+Double_t addX(Double_t x, THcShowerHit* h);
+Double_t addZ(Double_t x, THcShowerHit* h);
+Double_t addEpr(Double_t x, THcShowerHit* h);
+Double_t clX(THcShowerCluster* cluster);
+Double_t clZ(THcShowerCluster* cluster);
+Double_t clE(THcShowerCluster* cluster);
+Double_t clEpr(THcShowerCluster* cluster);
+Double_t clEplane(THcShowerCluster* cluster, Int_t iplane, Int_t side);
+
#endif
--
GitLab