From 871cd799458cd2e1013d30968b9094bf70501003 Mon Sep 17 00:00:00 2001
From: Vardan Tadevosyan <tadevosn@jlab.org>
Date: Wed, 12 Nov 2014 09:48:35 -0500
Subject: [PATCH] Add float precision calculation of the X coordinate of hit
 clusters in the HMS calorimeter.

    Add src/THcShowerCluster.h.float, a version of THcShowerCluster
    class with float precision calculation in clX method.
---
 src/THcShowerCluster.h.float | 254 +++++++++++++++++++++++++++++++++++
 1 file changed, 254 insertions(+)
 create mode 100644 src/THcShowerCluster.h.float

diff --git a/src/THcShowerCluster.h.float b/src/THcShowerCluster.h.float
new file mode 100644
index 0000000..796fb39
--- /dev/null
+++ b/src/THcShowerCluster.h.float
@@ -0,0 +1,254 @@
+#ifndef ROOT_THcShowerCluster
+#define ROOT_THcShowerCluster
+
+//HMS calorimeter hit cluster, version 3.
+
+//#include "THcShower.h"
+#include "THcShowerHit.h"
+
+//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.
+  //
+  Double_t clX() {
+    //    Double_t x_sum=0.;
+    //    Double_t Etot=0.;
+    float x_sum=0.;
+    float Etot=0.;
+
+    //    for (THcShowerHitIt it=THcShowerHitList::begin();
+    //	 it!=THcShowerHitList::end(); ++it) {
+
+    // Consistent with Engine.
+    //
+    for (THcShowerHitIt it=THcShowerHitList::end()-1;
+    	 it>=THcShowerHitList::begin(); it--) {
+      //x_sum += (*it)->hitX() * (*it)->hitE();
+      //Etot += (*it)->hitE();
+      x_sum += float((*it)->hitX()) * float((*it)->hitE());
+      Etot += float((*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();
+  }
+
+};
+
+//-----------------------------------------------------------------------------
+
+//Alias for container of clusters and for its iterator
+//
+typedef vector<THcShowerCluster*> THcShClusterList;
+typedef THcShClusterList::iterator THcShClusterIt;
+
+//List of clusters
+//
+class THcShowerClusterList : private THcShClusterList {
+
+ public:
+
+  THcShowerClusterList() {
+    //    cout << "Dummy THcShowerClusterList object created" << endl;
+  }
+
+  ~THcShowerClusterList() {
+    purge();
+  }
+
+  // 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();
+  }
+
+  //____________________________________________________________________________
+
+  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
+
+  }
+
+};
+
+#endif
-- 
GitLab