#include <math.h>
#include <iostream.h>
#include "dMutFitParams.h"
#include "dMutArrayParams.h"
#include "dMutCSCGeom.h"
#include "dMutGlobalPar.h"
#include "dMutClusterGhitRel.h"
#include "dMutClusterGhitRelN.h"
#include "mMutFitCathodeClustsT.h"
#include "dMutCathodeClustersN.h"
#include "dMutCathodeClusters.h"
#include "mumdigipar.h"
#include "PSTORE.h"
#include "sublink.h"
#include "fpmlink.h"
#include "MUMPAR.h"
#include "MutGeom.h"

#include "TNtuple.h"
extern TNtuple *MutClustDataNtuple;

//
//  Prototypes for functions called:
//
#include "mMutFitClusterP.h"
#include "mMutNestCathodeClustsCCP.h"
#include "mMutUnnestCathodeClustsCCP.h"
#include "mMutXYIntersectionP.h"

/**
  
   DESCRIPTION: Fit cathode clusters in digitization bank to 1 or more
     Mathieson distributions.  Choose the number to fit based on the
     width of the cathode cluster.  If poor chi-square result, try
     fitting a different number if there are enough strips in the 
     cluster.   <p>

   HISTORY:
               03dec98 Doug Fields (DEF) covnvert to C++
               00jan96-v000a-hpl- Created by stic Version   <p>
  
   @param   <p>
            IN:   <p>
             mumgeo     - Muon tracking geometry parameters from PISA   <p>
           mumgeo_h     -    <p>
         mumdigipar     - Input digitization parameters for tracking   <p>
       mumdigipar_h     -    <p>
         INOUT:   <p>
    dMutCathodeClusters - Table of found cathode clusters (sets of
                          contiguously hit strips)   <p>
    dMutCathodeClusters_h -    <p>
           OUT:   <p>

   @return     STAF Condition Value   <p>

   @author     Melynda Brooks, LANL, mbrooks@lanl.gov
*/
  
long mMutFitCathodeClustsT_(
			    TABLE_HEAD_ST          *mumgeo_h,        
			    MUMGEO_ST              *mumgeo,
			    TABLE_HEAD_ST          *mumdigipar_h,  
			    MUMDIGIPAR_ST          *mumdigipar ,
			    TABLE_HEAD_ST          *dMutCathodeClusters_h,
			    DMUTCATHODECLUSTERS_ST *dMutCathodeClusters,
			    TABLE_HEAD_ST          *dMutHistos_h,  
			    DMUTHISTOS_ST          *dMutHistos ,
			    TABLE_HEAD_ST          *dBbcOut_h,  
			    DBBCOUT_ST             *dBbcOut )
{
// 
//  Local variables:
//
  int goodnest;
  int stripfit;
  double a[MMAX];
  long int ioctant;
  long int ihalf;
  long int iplane;
  long int istation;
  long int iarm;
  long int igap;
  long int NumberOfGaps;
  long int NumberOfStrips;
  long int maxhit;  // true if fitted position is at edge of chamber
  long int icath;
  long int ihit;
  long int i;
  long int nclusts;
  long int iclust;
  long int cluswidf;
  long int cluswidfnew;
  long int istartf;
  int istripstart;
  const double SQRT12 = 3.464102;
  long int nhits;
  long int istrip;
  long int iflag;
  int fitcharge;			    	// Fit charge and centroid or just centroid
  double qmeas[MAXWID];
  double chisqnew=0; 
  double xbegin=0;
  double xbeginnew=0;
  double padwid=0;
  double sigq[MAXWID];         // Sigma of charge; large if overflow, etc.
  double xhit[MAXFIT];
  double qanode[MAXFIT];
  double covar[MMAX][MMAX];
  double xhitnew[MAXFIT];
  double qanodenew[MAXFIT];
  double chisqnewer=0;
  double qmaxnew=0;
  double xinit[MAXFIT];
  double qanodeinit[MAXFIT];
  double xtemp;
  double qtemp;
  double qmax;
  int found;
  long int temp;
  
  dMutCathodeClustersN_ST  *dMutCathodeClustersN[NumberOfArms]
                                                [NumberOfStations]
                                                [NumberOfOctants]
                                                [NumberOfHalfOctants]
                                                [3]
                                                [NumberOfCathodePlanes];
 
    MutArm* ArmPointer;
    int CathodePlane;
    double AnodeAngle, CathodeAngle, StereoAngle;

    float qpeak, qtot;
    float ntpvar[28];

#define PIBY2 1.570796
#define MAXWIDTH 25

// 
//  Other resolution degradation variables:
//

    static int EventNum = 0;
    int NumClus1, NumClus2, NumClus3;  //Number of clusters in each station
    int nclus;

//**************************************************************************************

    EventNum++;

//
//  Call routine to put cathode clusters into nested, multi-dimensional
//  structure, for easy looping.
//
    goodnest = mMutNestCathodeClustsCC( mumdigipar_h,
				       mumdigipar,
				       dMutCathodeClusters_h,
				       dMutCathodeClusters, 
				       dMutCathodeClustersN
				       );
//
//  Loop over each arm, station, and octant and fit cathode clusters to
//  1 or more Mathieson distributions. 
// 
      for (iarm = 0; iarm < NumberOfArms; ++iarm)
      { 
        switch (iarm){
          case (South):
            ArmPointer = SouthArm();
            break;
          case (North):
            ArmPointer = NorthArm();
            break;
          default:
            cout << "ERROR: invalid arm number used in mMutFitCathodeClustsT" << endl;
            ArmPointer = 0;
            break;
         }
        for (istation = 0; istation < NumberOfStations; ++istation)
        { 
          for (ioctant = 0; ioctant < NumberOfOctants; ++ioctant)
          { 
            for (ihalf= 0; ihalf< NumberOfHalfOctants; ++ihalf)
            { 
              NumberOfGaps = ArmPointer->f_pMutStations[istation]->
                             f_pMutOctants[ioctant]->
                             f_pMutHalfOctants[ihalf]->getNumberOfGaps();
               for (igap = 0; igap < NumberOfGaps; ++igap)
               { 
//
//  Look at both cathode planes
//
                for (icath = 0; icath < NumberOfCathodePlanes; ++icath)
                {
//
// Stereo angle for the plane
//
                   iplane = igap*2 + icath;
                   nclusts = dMutCathodeClustersN[iarm][istation][ioctant]
			       [ihalf][igap][icath]->nhits;

                   if (icath==0) CathodePlane = Cathode1;
                   else if (icath==1) CathodePlane = Cathode2;
                   else{
                     cout << "ERROR: icath not equal to 0 or 1 in mMutFitCathodeClustsT.cc" << endl;
                     CathodePlane = Cathode1;
                   }
                   AnodeAngle = ArmPointer->f_pMutStations[istation]->
                             f_pMutOctants[ioctant]->
                             f_pMutHalfOctants[ihalf]->
                             f_pMutGaps[igap]->
                             f_pMutPlanes[Wire]->
                             f_pMutWires[0]->getAngle();
                   CathodeAngle = ArmPointer->f_pMutStations[istation]->
                             f_pMutOctants[ioctant]->
                             f_pMutHalfOctants[ihalf]->
                             f_pMutGaps[igap]->
                             f_pMutPlanes[CathodePlane]->
                             f_pMutStrips[0]->getAngle();
                   StereoAngle = abs(AnodeAngle - CathodeAngle) - PIBY2;

//  Get number of strips in this plane so we can make sure that fitted positions
//  of clusters do not exceed the active area of the detector:

                   NumberOfStrips = ArmPointer->f_pMutStations[istation]->
                             f_pMutOctants[ioctant]->
                             f_pMutHalfOctants[ihalf]->
                             f_pMutGaps[igap]->
                             f_pMutPlanes[CathodePlane]->getNumElements() - 1;
//
//  Loop over all the clusters that have been found in the fine cathode
//  planes, and extract the cluster widths.  
//
                   for (iclust = 0; iclust < nclusts; ++iclust)
                   {   
                     cluswidf = dMutCathodeClustersN[iarm][istation][ioctant]
				  [ihalf][igap][icath]->clust[iclust]->cluswid;
                     istartf = dMutCathodeClustersN[iarm][istation][ioctant]
				  [ihalf][igap][icath]->clust[iclust]->stripbegin;
                     xbegin = (double)istartf;

                     if (cluswidf < MAXWIDTH) {
//
// Try to determine the number of tracks to fit:
//
                      nhits = 0;
                      qmax = 0;
                      for (i = 0; i < MAXFIT; ++i)
			          {
                         qanode[i] = 0.0;
                         xhit[i] = 0.0;
                      }

// Set some boundary conditions:

                      qmeas[0] = 0.0;
                      sigq[0] = 200.0;
                      qmeas[cluswidf+1] = 0.0;
                      sigq[cluswidf+1] = 200.0;

                      for (istrip = 0; istrip < cluswidf; ++istrip)
	              {
                         qmeas[istrip+1] = (double)dMutCathodeClustersN[iarm][istation]
			   [ioctant][ihalf][igap][icath]->clust[iclust]->qstrip[istrip];

                         sigq[istrip+1] = (double)dMutCathodeClustersN[iarm][istation]
                           [ioctant][ihalf][igap][icath]->clust[iclust]->qrms[istrip];

                         if (qmeas[istrip+1] > qmax) 
		         {
			       qmax = qmeas[istrip+1];
		         }
                         if (istrip == 0 && qmeas[istrip+1] > 
			         dMutCathodeClustersN[iarm][istation][ioctant]
				 [ihalf][igap][icath]->clust[iclust]->qstrip[istrip+1])
		         {
                           qanode[nhits] = qmeas[istrip+1];
                           xhit[nhits] = (double)(istartf + istrip);
                           nhits += 1;
		         }
                         else if (istrip == (cluswidf-1) && qmeas[istrip+1] > 
							 qmeas[istrip])
		         {
 		           found = 0;				  
                           for (i = 0; i < MAXFIT; ++i)
				           {
                              if (qanode[i] == 0 && found == 0)
					          {
                                found = 1;
                                qanode[nhits] = qmeas[istrip+1];
                                xhit[nhits] = (double)(istartf + istrip);
                                nhits += 1;
                              }
                           }
		         }
                         else if (qmeas[istrip+1] >= qmeas[istrip] && 
				 qmeas[istrip+1] > dMutCathodeClustersN[iarm]
				 [istation][ioctant][ihalf][igap][icath]->clust[iclust]->
				 qstrip[istrip+1])
                         {    
                           found = 0;
                           for (i = 0; i < MAXFIT; ++i)
			   {
                              if (qanode[i] == 0 && found == 0)
		              {
                                found = 1;
                                qanode[nhits] = qmeas[istrip+1];
                                xhit[nhits] = (double)(istartf + istrip);
                                nhits += 1;
			      }
			   }
			 }
		       }   // Loop over strips
//
// End of loop over strips in cluster
//
                      if (nhits > cluswidf/2) {
			  nhits = cluswidf/2;
		      }
                      if (nhits > MAXFIT) {
			  nhits = MAXFIT;
		      }
                      fitcharge = 1;
                      for (i = 0; i < nhits; ++i) {
                         xinit[i] = xhit[i];
                         qanodeinit[i] = qanode[i];
                      }
                      if (nhits > 0) {

                        cluswidfnew = cluswidf + 2;
                        xbeginnew = xbegin - 1;
                        padwid = ArmPointer->f_pMutStations[istation]->
                             f_pMutOctants[ioctant]->
                             f_pMutHalfOctants[ihalf]->
                             f_pMutGaps[igap]->
                             f_pMutPlanes[CathodePlane]->
                             f_pMutStrips[istartf]->getStripSpacing()/2.0;

                        for (i=0; i<nhits; i++){
                          a[i] =xhit[i];
                          a[nhits+i] = qanode[i];
                        }
                        istripstart = (int)istartf - 1;

                        stripfit = mMutFitCluster(qmeas, sigq, &(int)cluswidfnew,
                          &istripstart, &padwid, &(int)nhits, a, covar, &chisqnew);

                        for (i=0; i<nhits; i++){
                          xhit[i] = a[i];
                          qanode[i] = a[nhits+i];
                        }

//
//  If poor chi-square from fit, try alternative fits and see if
//  they give better chi-square result:
//
                        qmaxnew = 0.0;
                        if (chisqnew > 20) {
                          if (nhits < (cluswidf/2.0) && nhits < MAXFIT) {
                            for (i = 0; i < nhits; ++i) {
                               xhitnew[i] = xinit[i];
                               qanodenew[i] = qanodeinit[i]; 
		            }
                            qanodenew[nhits] = 0;
                            for (istrip = 0; istrip < cluswidf; ++istrip) {
                               found = 0;
                               for (i = 0; i < nhits; ++i) {
                                  if (xinit[i] == (istartf + istrip)) {
                                    found = 1;
                                  }
			       }
                               if (found == 0) {
                                 if (qmeas[istrip+1] > qanodenew[nhits]) {
                                   qanodenew[nhits] = qmeas[istrip+1];
                                   xhitnew[nhits] = istartf + istrip;
			         }
			       }
			    }  // Loop over strips
//
//  Order hits
//
                            ihit = 0;
                            while (ihit < nhits-1 && xhitnew[nhits] >
                                 xhitnew[ihit] && xhitnew[nhits] <
								 xhitnew[nhits-1])
			    {
			       ihit += 1;
		  	    }
                            if (ihit == 0 && xhitnew[nhits] < xhitnew[0]) {
			          ihit = -1;
		            }
                            else ihit=nhits-1;
                            xtemp = xhitnew[nhits];
                            qtemp = qanodenew[nhits];
                            for (i = nhits-1; i >= (ihit + 1); --i) {
                               xhitnew[i + 1] = xhitnew[i];
                               qanodenew[i + 1] = qanodenew[i];
                            }
                            xhitnew[ihit + 1] = xtemp;
                            qanodenew[ihit + 1] = qtemp;
//
// Do fitting
//
  			    temp = nhits+1;
                            cluswidfnew = cluswidf + 2;
                            xbeginnew = xbegin - 1;

                            for (i=0; i<nhits; i++){
                              a[i] = rint(xhitnew[i]);
                              a[nhits+i] = qanodenew[i];
                            }
                            istripstart = (int)istartf - 1;

                            stripfit = mMutFitCluster(qmeas, sigq, &(int)cluswidfnew,
                              &istripstart, &padwid, &(int)nhits, a, covar, &chisqnewer);


                            for (i=0; i<nhits; i++){
                              xhit[i] = a[i];
                              qanode[i] = a[nhits+i];
                            }

                            if (chisqnewer < chisqnew) {
                              nhits += 1;
                              chisqnew = chisqnewer;
                              for (ihit = 0; ihit < nhits; ++ihit)
					          {
                                 xhit[ihit] = xhitnew[ihit];
                                 qanode[ihit] = qanodenew[ihit];
                              }
		            }
                          }    // if fit fewer than could fit
                        }      // if poor chi-square from first fit

                        dMutCathodeClustersN[iarm][istation][ioctant]
                          [ihalf][igap][icath]->clust[iclust]->chisqr = chisqnew;
                        for (i = 0; i < nhits; ++i)
	                {

// If fitted position is beyond active area of chamber, truncate it to the edge:

                           maxhit = 0;
                           if (xhit[i] > NumberOfStrips){
                             xhit[i] = (float)NumberOfStrips;
                             maxhit = 1;
                           }
                           if (xhit[i] < 0){
                             xhit[i] = 0;
                             maxhit = 1;
                           }

                           dMutCathodeClustersN[iarm][istation][ioctant]
                             [ihalf][igap][icath]->clust[iclust]->xfit[i] 
                              = (float)xhit[i];
                           dMutCathodeClustersN[iarm][istation][ioctant]
                             [ihalf][igap][icath]->clust[iclust]->qfit[i] 
                             = (float)qanode[i];
                           if (nhits == 1)
		           {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = mumdigipar->resolution[istation]*
			         mumdigipar->wire_spacing[istation] + fabs(tan(StereoAngle)*
			         mumdigipar->wire_spacing[istation])*2.0/SQRT12;
			   }
                           else if (nhits == 2)
		           {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
		 	         clust[iclust]->res[i] = 1.1*mumdigipar->resolution[istation]*
			         mumdigipar->wire_spacing[istation] + fabs(tan(StereoAngle)*
			         mumdigipar->wire_spacing[istation])*2.0/SQRT12;
                           }
		           else
		           {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = 1.2*mumdigipar->resolution[istation]*
			         mumdigipar->wire_spacing[istation] + fabs(tan(StereoAngle)*
			         mumdigipar->wire_spacing[istation])*2.0/SQRT12;
		           }
//
//  Modify resolution if poor fit:
//
                           if (cluswidf/nhits < 2)
			   {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = 1.5*mumdigipar->resolution[istation]*
			         mumdigipar->wire_spacing[istation] + fabs(tan(StereoAngle)*
			         mumdigipar->wire_spacing[istation])*2.0/SQRT12;
			   }
                           if (chisqnew > 20) 
			   {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = cluswidf*mumdigipar->wire_spacing
				 [istation];
		           }
                           if (qmax > adcmax/adcconv - 1)
				           {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = cluswidf*mumdigipar->wire_spacing
				 [istation];
		           } 
                           if (maxhit) {
                             dMutCathodeClustersN[iarm][istation][ioctant]
                                 [ihalf][igap][icath]->
			         clust[iclust]->res[i] = cluswidf*mumdigipar->wire_spacing
				 [istation];
		           } 
	                }      
//
// End of loop over fitted positions
//
                        dMutCathodeClustersN[iarm][istation][ioctant]
                           [ihalf][igap][icath]->clust[iclust]->nfit = nhits;
	              }
                      else     
		      {

// Not enough info to fit cluster

  //  Make sure that stored hit position does not exceed the bounds of the chamber:

                        if (istartf + cluswidf/2.0 > NumberOfStrips){
                          dMutCathodeClustersN[iarm][istation][ioctant]
                             [ihalf][igap][icath]->clust[iclust]->xfit[0] 
                             = (float)NumberOfStrips;
                        }
                        else{
                          dMutCathodeClustersN[iarm][istation][ioctant]
                             [ihalf][igap][icath]->clust[iclust]->xfit[0] 
                             = istartf + cluswidf/2.0;
                        }
                        dMutCathodeClustersN[iarm][istation][ioctant]
                           [ihalf][igap][icath]->clust[iclust]->qfit[0] 
                           = qanode[0];
                        dMutCathodeClustersN[iarm][istation][ioctant]
                           [ihalf][igap][icath]->clust[iclust]->res[0] 
                           = cluswidf*mumdigipar->wire_spacing[istation];
                        dMutCathodeClustersN[iarm][istation][ioctant]
                           [ihalf][igap][icath]->clust[iclust]->chisqr = 100.;
                        dMutCathodeClustersN[iarm][istation][ioctant]
                           [ihalf][igap][icath]->clust[iclust]->nfit = 1;
		      }   // Test on cluster width
                    }  // if cluster width <= MAXWID
                    else{
                      cout << "Not fitting cluster which exceeds MAXWID" << endl;
                      dMutCathodeClustersN[iarm][istation][ioctant]
                         [ihalf][igap][icath]->clust[iclust]->xfit[0] 
                         = istartf + cluswidf/2.0;
                      dMutCathodeClustersN[iarm][istation][ioctant]
                         [ihalf][igap][icath]->clust[iclust]->qfit[0] 
                         = qanode[0];
                      dMutCathodeClustersN[iarm][istation][ioctant]
                         [ihalf][igap][icath]->clust[iclust]->res[0] 
                         = cluswidf*mumdigipar->wire_spacing[istation];
                      dMutCathodeClustersN[iarm][istation][ioctant]
                         [ihalf][igap][icath]->clust[iclust]->chisqr = 100.;
                      dMutCathodeClustersN[iarm][istation][ioctant]
                         [ihalf][igap][icath]->clust[iclust]->nfit = 1;
                    }
                  }      // Loop over fine cathode clusters
                }      // Loop over cathode planes
              }      // Loop over gaps
            }      // Loop over half octants
          }      // Loop over octants
        }      // Loop over stations
      }      // Loop over muon arms
 
//  Put nested structures back into linear structures, for now:

    goodnest = mMutUnnestCathodeClustsCC( mumdigipar_h,
				   mumdigipar,
				   dMutCathodeClusters_h,
                                   dMutCathodeClusters, 
				   dMutCathodeClustersN);

   if (dMutHistos->ClustData){
     NumClus1 = 0;
     NumClus2 = 0;
     NumClus3 = 0;
     nclus = 0;
     for (ihit=0; ihit<4; ihit++) ntpvar[20 + ihit] = -999;
     for (ihit=0; ihit<dMutCathodeClusters_h->nok; ihit++){
       if (dMutCathodeClusters[ihit].station==0) NumClus1++;
       if (dMutCathodeClusters[ihit].station==1) NumClus2++;
       if (dMutCathodeClusters[ihit].station==2) {
         NumClus3++;
         if (nclus<4){
           ntpvar[20+nclus] = dMutCathodeClusters[ihit].xfit[0];
           ntpvar[24+nclus] = dMutCathodeClusters[ihit].octant; 
           nclus++;
         }
       }
     }
     for (ihit=0; ihit<dMutCathodeClusters_h->nok; ihit++){
       qpeak = 0.0;
       qtot = 0.0;
       for (istrip=0; istrip<dMutCathodeClusters[ihit].cluswid; istrip++){
          qtot += dMutCathodeClusters[ihit].qstrip[istrip];
          if (dMutCathodeClusters[ihit].qstrip[istrip]>qpeak){
            qpeak = dMutCathodeClusters[ihit].qstrip[istrip];
          }
       }
       ntpvar[0] = (float)dMutCathodeClusters[ihit].arm;
       ntpvar[1] = (float)dMutCathodeClusters[ihit].station;
       ntpvar[2] = (float)dMutCathodeClusters[ihit].octant;
       ntpvar[3] = (float)dMutCathodeClusters[ihit].halfoctant;
       ntpvar[4] = (float)dMutCathodeClusters[ihit].gap;
       ntpvar[5] = (float)dMutCathodeClusters[ihit].cath;
       ntpvar[6] = (float)dMutCathodeClusters[ihit].stripbegin;
       ntpvar[7] = (float)dMutCathodeClusters[ihit].cluswid;
       ntpvar[8] = dMutCathodeClusters[ihit].xfit[0];
       ntpvar[9] = dMutCathodeClusters[ihit].qfit[0];
       ntpvar[10] = dMutCathodeClusters[ihit].chisqr; 
       ntpvar[11] = (float)dMutCathodeClusters[ihit].nfit;
       ntpvar[12] = (float)EventNum;
       ntpvar[13] = qpeak;
       ntpvar[14] = (float)NumClus1;
       ntpvar[15] = (float)NumClus2;
       ntpvar[16] = (float)NumClus3;
       ntpvar[17] = qtot;
       ntpvar[18] = (float)dBbcOut[0].NhitPmtSouth;
       ntpvar[19] = (float)dBbcOut[0].NhitPmtNorth;
       MutClustDataNtuple->Fill(ntpvar);
     }
   }  // If filling ClustData ntuple

// Delete nested structures:

    for (iarm=0; iarm<NumberOfArms; iarm++) {
      for (istation=0; istation<NumberOfStations; istation++) {
        for (ioctant=0; ioctant<NumberOfOctants; ioctant++) {
          for (ihalf=0; ihalf<NumberOfHalfOctants; ihalf++) {
            for (igap=0; igap<3; igap++) {
              for (icath=0; icath<NumberOfCathodePlanes; icath++) {
                if(dMutCathodeClustersN[iarm][istation][ioctant]
                   [ihalf][igap][icath]){
                  for (ihit=0; ihit<MUT_MAX_CLHITS; ihit++){
                    if(dMutCathodeClustersN[iarm][istation][ioctant]
                       [ihalf][igap][icath]->clust[ihit]){
                      delete dMutCathodeClustersN[iarm][istation][ioctant]
                        [ihalf][igap][icath]->clust[ihit];
                    }
                  }
                  delete dMutCathodeClustersN[iarm][istation][ioctant]
                    [ihalf][igap][icath];
                  dMutCathodeClustersN[iarm][istation][ioctant]
                    [ihalf][igap][icath] = NULL;
                }  //If pointer is non-zero
              }  //Loop over planes
            }  //Loop over gaps
          }  //Loop over half octants
        }  //Loop over octants
      }  //Loop over stations
    }  //Loop over arms


return STAFCV_OK;
}