diff --git a/src/THcShower.cxx b/src/THcShower.cxx
index a001241ecb04ae6c8f78beb1c6d4a4730bc56093..6de4c5cfdbc03d19e24cd2a64b7f30ae247091aa 100644
--- a/src/THcShower.cxx
+++ b/src/THcShower.cxx
@@ -40,6 +40,8 @@ THcShower::THcShower( const char* name, const char* description,
{
// Constructor
fNLayers = 0; // No layers until we make them
+
+ fClusterList = new THcShowerClusterList;
}
//_____________________________________________________________________________
@@ -507,6 +509,8 @@ void THcShower::Clear(Option_t* opt)
fTREpl_neg_cor[ip] = -0.;
}
+ fClusterList->clear();
+
}
//_____________________________________________________________________________
@@ -645,10 +649,10 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
// Create list of clusters and fill it.
- THcShowerClusterList* ClusterList = new THcShowerClusterList;
- ClusterList->ClusterHits(HitList);
+ // fClusterList = new THcShowerClusterList; //shall be allocated before
+ fClusterList->ClusterHits(HitList);
- fNclust = (*ClusterList).NbClusters(); //number of clusters
+ fNclust = (*fClusterList).NbClusters(); //number of clusters
//Print out the cluster list.
@@ -658,7 +662,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
for (Int_t i=0; i!=fNclust; i++) {
- THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+ THcShowerCluster* cluster = (*fClusterList).ListedCluster(i);
//cout << "Cluster #" << i
//<<": E=" << (*cluster).clE()
@@ -683,13 +687,13 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
if (fNclust != 0) {
- THcShowerCluster* MaxCluster = (*ClusterList).ListedCluster(0);
+ THcShowerCluster* MaxCluster = (*fClusterList).ListedCluster(0);
fE = (*MaxCluster).clE();
for (Int_t i=1; i<fNclust; i++) {
- THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+ THcShowerCluster* cluster = (*fClusterList).ListedCluster(i);
Double_t E = (*cluster).clE();
@@ -738,7 +742,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
// Associate a cluster to the track.
- Int_t mclust = MatchCluster(theTrack, ClusterList, Xtr, Ytr);
+ Int_t mclust = MatchCluster(theTrack, Xtr, Ytr);
// if (mclust >= 0) fNtracks++; // number of shower tracks (This is not consistent with engine!)
@@ -752,7 +756,7 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
if (mclust >= 0) { // if there is a matched cluster
//Assign energies (not Y corrected) of the matched cluster to the track.
- THcShowerCluster* cluster = (*ClusterList).ListedCluster(mclust);
+ THcShowerCluster* cluster = (*fClusterList).ListedCluster(mclust);
fTRE = (*cluster).clE();
fTREpr = (*cluster).clEpr();
@@ -820,7 +824,6 @@ Int_t THcShower::CoarseProcess( TClonesArray& tracks)
//-----------------------------------------------------------------------------
Int_t THcShower::MatchCluster(THaTrack* Track,
- THcShowerClusterList* ClusterList,
Double_t& XTrFront, Double_t& YTrFront)
{
// Match a cluster to a given track. Return the cluster number,
@@ -902,7 +905,7 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
for (Int_t i=0; i<fNclust; i++) {
- THcShowerCluster* cluster = (*ClusterList).ListedCluster(i);
+ THcShowerCluster* cluster = (*fClusterList).ListedCluster(i);
Double_t dx = TMath::Abs( (*cluster).clX() - XTrFront );
@@ -923,6 +926,72 @@ Int_t THcShower::MatchCluster(THaTrack* Track,
return mclust;
}
+//_____________________________________________________________________________
+Float_t THcShower::GetShEnergy(THaTrack* Track) {
+
+ // Get total energy deposited in the cluster matched to the given
+ // spectrometer Track.
+
+ // Track coordinates at front of the calorimeter, initialize out of
+ // acceptance.
+ Double_t Xtr = -75.;
+ Double_t Ytr = -40.;
+
+ // Associate a cluster to the track.
+
+ Int_t mclust = MatchCluster(Track, Xtr, Ytr);
+
+ if (fdbg_tracks_cal) {
+ cout << "GetShEnergy: Track X = " << Xtr << " Y = " << Ytr;
+ if (mclust >= 0)
+ cout << " matched cluster # " << mclust << endl;
+ else
+ cout << " no matched cluster found" << endl;
+ }
+
+ // Get matched cluster.
+ THcShowerCluster* cluster = (*fClusterList).ListedCluster(mclust);
+
+ // Coordinate corrected total energy deposition in the cluster.
+
+ Float_t Etrk = 0.;
+ if (mclust >= 0) { // if there is a matched cluster
+
+ // Correct track energy depositions for the impact coordinate.
+
+ for (Int_t ip=0; ip<fNLayers; ip++) {
+
+ // Coordinate correction factors for positive and negative sides,
+ // different for single PMT counters in the 1-st two layes and for
+ // 2 PMT counters in the rear two layers.
+ Float_t corpos = 1.;
+ Float_t corneg = 1.;
+ if (ip < fNegCols) {
+ corpos = Ycor(Ytr,0);
+ corneg = Ycor(Ytr,1);
+ }
+ else {
+ corpos = Ycor(Ytr);
+ corneg = 0.;
+ }
+
+ if (fdbg_tracks_cal) {
+ cout << " Plane " << ip << " Ytr = " << Ytr
+ << " corpos = " << corpos
+ << " corneg = " << corneg << endl;
+ }
+
+ Etrk += (*cluster).clEplane(ip,0) * corpos;
+ Etrk += (*cluster).clEplane(ip,1) * corneg;;
+
+ } //over planes
+
+ } //mclust>=0
+
+ if (fdbg_tracks_cal) cout << " Etrk = " << Etrk << endl;
+
+ return Etrk;
+}
//_____________________________________________________________________________
Int_t THcShower::FineProcess( TClonesArray& tracks )
diff --git a/src/THcShower.h b/src/THcShower.h
index 04b2a3a2f6b069751e8f4dfd4c157f3a1237cfa5..04d351dfd28b55b8969390a3188663f55b88dcdc 100644
--- a/src/THcShower.h
+++ b/src/THcShower.h
@@ -109,6 +109,11 @@ public:
return (fCcor + sign*y)/(fCcor + sign*y/fDcor);
}
+ // Get total energy deposited in the cluster matched to the given
+ // spectrometer Track.
+
+ Float_t GetShEnergy(THaTrack*);
+
THcShower(); // for ROOT I/O
protected:
@@ -137,6 +142,8 @@ protected:
// Int_t fNblk; // Number of blocks in main cluster
// Double_t* fEblk; // Energies of blocks in main cluster
+ THcShowerClusterList* fClusterList; // List of hit clusters
+
// Track quantities, are here for test purposes. Better to move
// to the track class(es) later on.
@@ -195,7 +202,7 @@ protected:
// cluster to track association method.
- Int_t MatchCluster(THaTrack*, THcShowerClusterList*, Double_t&, Double_t&);
+ Int_t MatchCluster(THaTrack*, Double_t&, Double_t&);
friend class THcShowerPlane; //to access debug flags.
diff --git a/src/THcShowerCluster.h b/src/THcShowerCluster.h
index f716e73c81d964ed98d7055ace40d2feaae076c4..59404365e0c74a6f4a3699112da5c311c3014428 100644
--- a/src/THcShowerCluster.h
+++ b/src/THcShowerCluster.h
@@ -169,6 +169,12 @@ class THcShowerClusterList : private THcShClusterList {
}
}
+ // Clear cluster list
+ //
+ void clear() {
+ THcShClusterList::clear();
+ }
+
//Put a cluster in the cluster list
//
void grow(THcShowerCluster* cluster) {