diff -N ClusterShapes.cpp
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ClusterShapes.cpp	21 May 2010 21:12:48 -0000	1.1
@@ -0,0 +1,118 @@

+/* -*- Mode: C++; tab-width: 2; indent-tabs-mode: nil; c-basic-offset: 2 -*- */
+
+#include "ClusterShapes.h"
+
+
+
+// #################################################
+// #####                                       #####
+// #####  Additional Structures and Functions  #####
+// #####                                       #####
+// #################################################
+//=============================================================================
+
+struct data {
+  int n;
+  float* x;
+  float* y;
+  float* z;
+  float* ex;
+  float* ey;
+  float* ez;
+};
+
+
+// ##########################################
+// #####                                #####
+// #####   Constructor and Destructor   #####
+// #####                                #####
+// ##########################################
+
+//=============================================================================
+
+ClusterShapes::ClusterShapes(int nhits, float* a, float* x, float* y, float* z){
+
+  _nHits = nhits;
+  _aHit = new float[_nHits];
+  _xHit = new float[_nHits];
+  _yHit = new float[_nHits];
+  _zHit = new float[_nHits];
+  _exHit = new float[_nHits];   
+  _eyHit = new float[_nHits];
+  _ezHit = new float[_nHits];         
+  _types = new int[_nHits];
+
+  _xl = new float[_nHits];
+  _xt = new float[_nHits];
+
+  _t = new float[_nHits];
+  _s = new float[_nHits]; 
+
+  for (int i(0); i < nhits; ++i) {
+    _aHit[i] = a[i];
+    _xHit[i] = x[i];
+    _yHit[i] = y[i];
+    _zHit[i] = z[i];
+    _exHit[i] = 1.0;
+    _eyHit[i] = 1.0;
+    _ezHit[i] = 1.0;
+    _types[i] = 1; // all hits are assumed to be "cyllindrical"
+  }
+
+  _ifNotGravity     = 1;
+  _ifNotWidth       = 1;
+  _ifNotInertia     = 1;
+  _ifNotEigensystem = 1;
+
+}
+
+
+//=============================================================================
+
+ClusterShapes::~ClusterShapes() {
+
+  delete[] _aHit;
+  delete[] _xHit;
+  delete[] _yHit;
+  delete[] _zHit;
+  delete[] _exHit;
+  delete[] _eyHit;
+  delete[] _ezHit;
+  delete[] _types;
+  delete[] _xl;
+  delete[] _xt;
+  delete[] _t;
+  delete[] _s;
+}
+
+//=============================================================================
+
+float* ClusterShapes::getCentreOfGravity() {
+  if (_ifNotGravity == 1) findGravity() ;
+  return &_analogGravity[0] ;
+}
+
+//=============================================================================
+
+void ClusterShapes::findGravity() 
+{
+    _totAmpl = 0. ;
+    for (int i(0); i < 3; ++i) 
+    {
+      _analogGravity[i] = 0.0 ;
+    }
+    for (int i(0); i < _nHits; ++i) {
+      _totAmpl+=_aHit[i] ;
+      _analogGravity[0]+=_aHit[i]*_xHit[i] ;
+      _analogGravity[1]+=_aHit[i]*_yHit[i] ;
+      _analogGravity[2]+=_aHit[i]*_zHit[i] ;
+    }
+    for (int i(0); i < 3; ++i) {
+      _analogGravity[i]/=_totAmpl ;
+    }
+    _xgr = _analogGravity[0];
+    _ygr = _analogGravity[1];
+    _zgr = _analogGravity[2];
+    _ifNotGravity = 0;
+}
+

diff -N ClusterShapes.h
--- /dev/null	1 Jan 1970 00:00:00 -0000
+++ ClusterShapes.h	21 May 2010 21:12:48 -0000	1.1
@@ -0,0 +1,113 @@

+#ifndef ClusterShapes_h
+#define ClusterShapes_h
+
+#include <stdio.h>
+#include <iostream>
+#include <iomanip>
+#include <string>
+#include <sstream>
+#include <cstdlib>
+
+#include <math.h>
+
+/**
+ *    Utility class to derive properties of clusters, such as centre of gravity,
+ *    axes of inertia, fits of the cluster shape and so on. All the details are
+ *    explained in the documentation of the methods. Several classes of the GSL 
+ *    (GNU Scientific Library) are needed in this class.
+ *
+ *    @authors V. Morgunov (ITEP/DESY), A. Raspereza (DESY), O. Wendt (DESY)
+ *    @version $Id: ClusterShapes.h,v 1.1 2010/05/21 21:12:48 jeremy Exp $
+ *
+ *    --
+ *
+ *    Ported from MarlinUtil.  Only took findGravity() / getCentreOfGravity() code.
+ *    No GSL dependency.  --JM
+ *
+ */
+class ClusterShapes 
+{
+
+public:
+
+    /**
+     *    Constructor
+     *    @param nhits : number of hits in the cluster
+     *    @param a     : amplitudes of elements ('cells') of the cluster. Stored in 
+     *                   an array, with one entry for each element ('cell'). Each entry 
+     *                   is depending on coordinates x,y,z (Cartesian), which are stored 
+     *                   in the arrays x,y,z.
+     *    @param x,y,z : array of coordinates corresponding to the array of amplitudes a.
+     *
+     *
+     */
+    ClusterShapes(int nhits, float* a, float* x, float* y, float* z);
+
+
+    /**
+     *    Destructor
+     */
+    ~ClusterShapes();
+
+    /**
+     * returns an array, which represents a vector from the origin of the
+     * coordiante system, i.\ e.\ IP, to the centre of gravity of the cluster. The centre 
+     * of gravity is calculated with the energy of the entries of the cluster.
+     */
+    float* getCentreOfGravity();
+
+    /** US spelling of getCentreOfGravity */
+    inline float* getCenterOfGravity() { return getCentreOfGravity() ; }
+
+private:
+
+    void  findGravity();
+
+private:
+
+    int _nHits;
+
+    float* _aHit;
+    float* _xHit;
+    float* _yHit;
+    float* _zHit;
+    float* _exHit;
+    float* _eyHit;
+    float* _ezHit;
+    int* _types;
+    float* _xl;
+    float* _xt;
+    float* _t;
+    float* _s;
+
+    int   _ifNotGravity;
+    float _totAmpl;
+    float _radius;
+    float _xgr;
+    float _ygr;
+    float _zgr;
+    float _analogGravity[3];
+
+    int   _ifNotWidth;
+    float _analogWidth;
+
+    int   _ifNotInertia;
+    float _ValAnalogInertia[3];
+    float _VecAnalogInertia[9];
+
+    int _ifNotEigensystem;
+
+    int   _ifNotElipsoid;
+    float _r1           ;  // Cluster spatial axis length -- the largest
+    float _r2           ;  // Cluster spatial axis length -- less
+    float _r3           ;  // Cluster spatial axis length -- less
+    float _vol          ;  // Cluster ellipsoid volume
+    float _r_ave        ;  // Cluster average radius  (qubic root)
+    float _density      ;  // Cluster density
+    float _eccentricity ;  // Cluster Eccentricity
+    float _r1_forw      ;
+    float _r1_back      ; 
+
+};
+
+#endif