---

SiClusters.cxx

Go to the documentation of this file.

// $Id: SiClusters.cxx,v 1.3 2000/12/11 16:38:45 burnett Exp $

#include "reconstruction/SiClusters.h"
#include "DetectorData.h"

#include <float.h> // for FLT_MAX

using namespace std;

void SiClusters::Tray::flagHit (unsigned int idx)
{
    unsigned nx = xhitList.size();
    unsigned ix = idx%nx;
    unsigned iy = idx/nx;

    xhitList[ix]->flag = 1;
    yhitList[iy]->flag = 1;
}

void SiClusters::Tray::flagHit (enum SiData::Axis , unsigned int ) {}

void SiClusters::Tray::unflagHit (unsigned int idx)
{
    unsigned nx = xhitList.size();
    unsigned ix = idx%nx;
    unsigned iy = idx/nx;

    xhitList[ix]->flag = 0;
    yhitList[iy]->flag = 0;
}

void SiClusters::Tray::unflagHit (enum SiData::Axis , unsigned int ) {}

void SiClusters::Tray::clearHitFlags ()
{
     unsigned i;
     for( i =0; i<xhitList.size(); i++)
         xhitList[i]->flag = 0;
     for( i =0; i<yhitList.size(); i++)
         yhitList[i]->flag = 0;
}


SiClusters::SiClusters (int ntrays)
{
    trayList.resize(ntrays);
}

SiClusters::~SiClusters ()
{
    clear();
}

void SiClusters::loadFrom (const SiData& data)
{
    for(unsigned i=0; i<trayList.size(); i++) {

  //  First do x....
        int nx = data.nHits(SiData::X, i);
        if(nx) {

   // Find the start and end of a modules worth of data and  transfer
           unsigned lastMod = 9999999;
           int start=0, stop;
           for(stop=0; stop<nx; stop++) {
              if(lastMod != data.moduleId(SiData::X, i, stop)) {
                if((stop-start) > 0)
                  transferData(SiData::X, data, i, start, stop);
                lastMod = data.moduleId(SiData::X, i, stop);
                start = stop;
              }
            }
            if((stop-start) > 0) {
               transferData(SiData::X, data, i, start, stop);
            }
         }

  // Now do y....
        int ny = data.nHits(SiData::Y, i);
        if(ny) {
           unsigned lastMod = 9999999;
           int start=0, stop;
           for(stop=0; stop<ny; stop++) {
              if(lastMod != data.moduleId(SiData::Y, i, stop)) {
                if((stop-start) > 0)
                  transferData(SiData::Y, data, i, start, stop);
                lastMod = data.moduleId(SiData::Y, i, stop);
                start = stop;
              }
            }
            if((stop-start) > 0) {
               transferData(SiData::Y, data, i, start, stop);
            }
         }
      }

}

void SiClusters::clear ()
{
  for(unsigned i=0; i<trayList.size(); i++) {
    vector<Tray *> &moduleTrays = trayList[i];
    vector<Tray *>::iterator it = moduleTrays.begin();
    for ( ; it != moduleTrays.end(); it++ ) {
      delete *it;
    }
    moduleTrays.clear();
  }
}

int SiClusters::nHits (enum SiData::Axis a, int trayNum) const
{
  int count = 0;
  if ((trayNum < 0) || (trayNum >= trayList.size())) return 0;
  const vector<Tray *> &moduleTrays = trayList[trayNum];
  vector<Tray *>::const_iterator it = moduleTrays.begin();
  for( ; it != moduleTrays.end(); ++it ) {
    Tray *tray = *it;
      if(a==SiData::X) count += tray->xhitList.size();
      else             count += tray->yhitList.size();
  }
  return count;
}

Point SiClusters::nextHit (int trayN, int* index)
{
     static int trayNum = 999;
     static int modNum  = 999;

     static int ix = 0, iy = 0, nx = 0, ny = 0, pairs = 0;
     static Tray *tray = 0;

     *index = -1;           // return of negative signals no hits available
     if(trayN != trayNum) { // Note: trayN < 0 forces initialization
        modNum = 0;
        ix = iy = nx = ny = pairs = 0;
        if(trayN < 0)  {
           trayNum = 999;
           return Point(0., 0., 0.);
        }
        else
           trayNum = trayN;
     }

     if(pairs >= nx*ny) {
       // exhausted all combinations: next Module or terminate
        nx = ny = 0;
        vector <Tray *> &moduleTrays = trayList[trayNum];
        while(!nx || !ny) {
           if(modNum >= (int)moduleTrays.size()) return Point(0., 0., 0.);
           tray = moduleTrays[modNum++];
           nx = tray->xhitList.size();
           ny = tray->yhitList.size();
        }
        ix = iy = pairs = 0;
     }
     Point xpt , ypt;
     if(!tray->xhitList[ix]->flag)
        xpt = tray->xhitList[ix]->pos;
     else {
        pairs++;
        ix++;
        if(ix >= nx) {
           iy++;
           ix = 0;
           return nextHit(trayNum, index);
        }
     }

     if(!tray->yhitList[iy]->flag)
        ypt = tray->yhitList[iy]->pos;
     else {
        pairs += nx;
        iy++;
        ix = 0;
        return nextHit(trayNum, index);
     }

     float x = xpt.x();
     float y = ypt.y();
     float z = (xpt.z() + ypt.z())/2.;

     *index = pairs + 10000*trayNum + 1000000*tray->module;

     pairs++;
     ix++;
     if(ix >= nx) {
        iy++;
        ix = 0;
     }

     return Point(x, y, z);
}

Point SiClusters::nextHit (enum SiData::Axis a, int trayN, int* index)
{
     static int trayNum = 999;
     static int modNum  = 999;

     static unsigned int ix = 0,  nx = 0;
     static Tray *tray = 0;
     static vector<Cluster *> *hitList= 0;

     *index = -1;           // return of negative signals no hits available
     if(trayN != trayNum) { // Note: trayN < 0 forces initialization
        modNum = 0;
        ix  = nx = 0;
        if(trayN < 0)  {
           trayNum = 999;
           return Point(0., 0., 0.);
        }
        else
           trayNum = trayN;
     }

     if(ix >= nx) {   // exhausted all hits: next Module or terminate
        nx  = 0;
        vector<Tray *> &moduleTrays = trayList[trayNum];
        while(!nx) {
           if(modNum >= (int)moduleTrays.size()) return Point(0., 0., 0.);
           tray = moduleTrays[modNum++];

           if(a==SiData::X) hitList = &tray->xhitList;
           else             hitList = &tray->yhitList;
           nx = hitList->size();
        }
        ix = 0;
     }
     if ( ix >= hitList->size() ) {
       // this should never happen, but it did under Linux with protons
       *index = -1;
       return Point(0., 0., 0.);
     }
     Point xpt;
     if(!hitList->operator[](ix)->flag)  xpt = hitList->operator[](ix)->pos;
     else {
           ix++;
           return nextHit(a, trayNum, index);
     }

     *index = ix + 10000*trayNum + 1000000*tray->module;

     ix++;

     return xpt;
}

Point SiClusters::nearestHitInside (unsigned int trayNum, float dR,
                                    const Point& x0, int* index)
{
   float  x = 0., y = 0., z = 0., test= fabs(dR);

    *index = -1;
    if(test < 100000.) test *= test;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
          float xtmp2 = x0.x() - tray->xhitList[ix]->pos.x();
          xtmp2 *=  xtmp2;
          if(tray->xhitList[ix]->flag || xtmp2 > test) continue;

          for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
             float ytmp2 = x0.y() - tray->yhitList[iy]->pos.y();
             ytmp2 *= ytmp2;
             if(tray->yhitList[iy]->flag || ytmp2 > test) continue;

             if((xtmp2 + ytmp2) < test) {
                test = (xtmp2 + ytmp2);
                x = tray->xhitList[ix]->pos.x();
                y = tray->yhitList[iy]->pos.y();
                z = (tray->yhitList[iy]->pos.z()
                   + tray->xhitList[ix]->pos.z())/2.;
                *index = 10000*trayNum + 1000000*tray->module;
                *index += ix + iy*tray->xhitList.size();
             }
          }
       }
    }
    return Point(x, y, z);
}

Point SiClusters::nearestHitInside (unsigned int trayNum, float dX, float dY,
                                    const Point& x0, int* index)
{
   float  x = 0., y = 0., z = 0., test= FLT_MAX;

    *index = -1;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
          float xtmp = x0.x() - tray->xhitList[ix]->pos.x();
          if(tray->xhitList[ix]->flag || xtmp > dX) continue;

          for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
             float ytmp = x0.y() - tray->yhitList[iy]->pos.y();
             if(tray->yhitList[iy]->flag || ytmp > dY) continue;

             float dsq = xtmp*xtmp + ytmp*ytmp;
             if(dsq < test) {
                test = dsq;
                x = tray->xhitList[ix]->pos.x();
                y = tray->yhitList[iy]->pos.y();
                z = (tray->yhitList[iy]->pos.z()
                   + tray->xhitList[ix]->pos.z())/2.;
                *index = 10000*trayNum + 1000000*tray->module;
                *index += ix + iy*tray->xhitList.size();
              }
          }
       }
    }
    return Point(x, y, z);
}

Point SiClusters::nearestHitInside (enum SiData::Axis a, unsigned int trayNum,
                                    float dX, const Point& x0, int* index)
{
   static  vector<Cluster *> *hitList = 0;
   float  test= FLT_MAX;
   Point hit(0.,0.,0.);

    *index = -1;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       if(a==SiData::X) hitList = &tray->xhitList;
       else             hitList = &tray->yhitList;

       for(unsigned ix=0; ix<hitList->size(); ix++) {
          float xtmp;
          if(a==SiData::X) xtmp = fabs(x0.x()
                                  - hitList->operator[](ix)->pos.x());
          else             xtmp = fabs(x0.y()
                                  - hitList->operator[](ix)->pos.y());
  //      if(hitList->operator[](ix)->flag || xtmp > dX) continue;
          if( xtmp > dX) continue;

          if(xtmp < test) {
             test = xtmp;
             hit = hitList->operator[](ix)->pos;
             *index = 10000*trayNum + 1000000*tray->module;
             *index += ix;
          }
       }
    }
    return hit;
}

Point SiClusters::nearestHitOutside (unsigned int trayNum, float dR,
                                     const Point& x0, int* index)
{
   float  x = 0., y = 0., z = 0., test= FLT_MAX, dR2 = dR*dR;

   if(dR < 0) dR2 = -dR2;
    *index = -1;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
          if(tray->xhitList[ix]->flag) continue;
          float xtmp2 = x0.x() - tray->xhitList[ix]->pos.x();
          xtmp2 *=  xtmp2;

          for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
             if(tray->yhitList[iy]->flag) continue;
             float ytmp2 = x0.y() - tray->yhitList[iy]->pos.y();
             ytmp2 *= ytmp2;

             if((xtmp2 + ytmp2) < test && (xtmp2 + ytmp2) > dR2) {
                test = (xtmp2 + ytmp2);
                x = tray->xhitList[ix]->pos.x();
                y = tray->yhitList[iy]->pos.y();
                z = (tray->yhitList[iy]->pos.z()
                   + tray->xhitList[ix]->pos.z())/2.;
                *index = 10000*trayNum + 1000000*tray->module;
                *index += ix + iy*tray->xhitList.size();
             }
          }
       }
    }
    return Point(x, y, z);
}

Point SiClusters::nearestHitOutside(enum SiData::Axis a, unsigned int trayNum,
                                    float dX, const Point& x0, int* index)
{
   static  vector<Cluster *> *hitList = 0;
   float  test= FLT_MAX;
   Point hit(0.,0.,0.);

    *index = -1;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       if(a==SiData::X) hitList = &tray->xhitList;
       else             hitList = &tray->yhitList;

       for(unsigned ix=0; ix<hitList->size(); ix++) {
          if(hitList->operator[](ix)->flag) continue;
          float xtmp;
          if(a==SiData::X) xtmp = fabs(x0.x()
                                - hitList->operator[](ix)->pos.x());
          else             xtmp = fabs(x0.y()
                                - hitList->operator[](ix)->pos.y());

          if(xtmp < test && xtmp > dX) {
             test = xtmp;
             hit = hitList->operator[](ix)->pos;
             *index = 10000*trayNum + 1000000*tray->module;
             *index += ix;
          }
       }
    }
    return hit;
}

Point SiClusters::meanHit (unsigned int trayNum, float* rmsX, float* rmsY)
{
  double x = 0., y = 0., z = 0., x2 = 0., y2 = 0.,
        xnum = 0., ynum = 0., znum = 0.;

    *rmsX = 0.;
    *rmsY = 0.;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
          x += tray->xhitList[ix]->pos.x();
          z += tray->xhitList[ix]->pos.z();
          x2 += tray->xhitList[ix]->pos.x()*tray->xhitList[ix]->pos.x();
          xnum += 1.;
          znum += 1.;
       }

       for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
          y += tray->yhitList[iy]->pos.y();
          z += tray->yhitList[iy]->pos.z();
          y2 += tray->yhitList[iy]->pos.y()*tray->yhitList[iy]->pos.y();
          ynum += 1.;
          znum += 1.;
       }
    }
    if(xnum > 1) {
        x /= xnum;
        x2 /= xnum;
        double t = x2-x*x;
        *rmsX = t>0? sqrt(t):0;
    }
    if(ynum > 1) {
        y /= ynum;
        y2 /= ynum;
        double t = y2-y*y;
        *rmsY = t>0? sqrt(t):0;
    }
    if(xnum+ynum > 0.) z /= (xnum+ynum);

    return Point(x, y, z);
}
//########################################################################
Point SiClusters::meanHit (SiData::Axis axis, unsigned int trayNum, float* rms)
//########################################################################
{
  double x = 0., y = 0., z = 0., x2 = 0., y2 = 0.,
        xnum = 0., ynum = 0., znum = 0.;

    *rms = 0.;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       if (axis == SiData::X) {
           for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
               x += tray->xhitList[ix]->pos.x();
               z += tray->xhitList[ix]->pos.z();
               x2 += tray->xhitList[ix]->pos.x()*tray->xhitList[ix]->pos.x();
               xnum += 1.;
               znum += 1.;
           }
       }

       if (axis == SiData::Y) {
           for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
               y += tray->yhitList[iy]->pos.y();
               z += tray->yhitList[iy]->pos.z();
               y2 += tray->yhitList[iy]->pos.y()*tray->yhitList[iy]->pos.y();
               ynum += 1.;
               znum += 1.;
           }
       }
    }
    if(axis == SiData::X && xnum > 0) {
        x /= xnum;
        x2 /= xnum;
        double t = x2 - x*x;
        *rms = t>0? sqrt(t):0;
    }
    if (axis == SiData::Y && ynum > 0) {
        y /= ynum;
        y2 /= ynum;
        double t = y2 - y*y;
        *rms = t>0? sqrt(t):0;
    }
    if(xnum+ynum > 0.) z /= (xnum+ynum);

    return Point(x, y, z);
}
Point SiClusters::meanHitInside (unsigned int trayNum, float dXY,
                                 const Point& x0)
{
   float  x = 0., y = 0., z = 0.,  num = 0.;

   float dXYsq = dXY*dXY;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

       for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
          float xtmp2 = x0.x() - tray->xhitList[ix]->pos.x();
          xtmp2 *=  xtmp2;
          if(tray->xhitList[ix]->flag || xtmp2 > dXYsq) continue;

          for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
             float ytmp2 = x0.y() - tray->yhitList[iy]->pos.y();
             ytmp2 *= ytmp2;
             if(tray->yhitList[iy]->flag || ytmp2 > dXYsq) continue;

             if((xtmp2 + ytmp2) < dXYsq) {
                x += tray->xhitList[ix]->pos.x();
                y += tray->yhitList[iy]->pos.y();
                z += (tray->yhitList[iy]->pos.z()
                    + tray->xhitList[ix]->pos.z())/2.;
                num += 1.;
             }
          }
       }
    }
    if(num > 0) return Point(x/num, y/num, z/num);
    else return Point(0., 0., 0.);
}

int SiClusters::numberOfHitsNear (unsigned int trayNum, float dXY,
                                  const Point& x0)
{
    return numberOfHitsNear(trayNum, dXY, dXY, x0);
}

int SiClusters::numberOfHitsNear (unsigned int trayNum, float dX, float dY,
                                  const Point& x0)
{
  int  xnum = 0, ynum = 0;
  float d2_width = DetectorData::Tower::mod_width/2.;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

        for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
           if(fabs(x0.x() - tray->xhitList[ix]->pos.x()) < dX &&
              fabs(x0.y() - tray->xhitList[ix]->pos.y()) < d2_width)  xnum++;
        }

        for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
           if(fabs(x0.y() - tray->yhitList[iy]->pos.y()) < dY &&
              fabs(x0.x() - tray->yhitList[iy]->pos.x()) < d2_width)  ynum++;
       }
    }
    return xnum+ynum;
}
int SiClusters::numberOfHitsNear (enum SiData::Axis axis, unsigned int trayNum, float dR, const Point& x0)
{
  int  xnum = 0, ynum = 0;
  float d2_width = DetectorData::Tower::mod_width/2.;

    for(unsigned j=0; j<trayList[trayNum].size(); j++) {
       Tray *tray = trayList[trayNum][j];

                if (axis == SiData::X) {
                        for(unsigned ix=0; ix<tray->xhitList.size(); ix++) {
                        if(fabs(x0.x() - tray->xhitList[ix]->pos.x()) < dR &&
                                  fabs(x0.y() - tray->xhitList[ix]->pos.y()) < d2_width)  xnum++;
                        }
                } else {
                        for(unsigned iy=0; iy<tray->yhitList.size(); iy++) {
                        if(fabs(x0.y() - tray->yhitList[iy]->pos.y()) < dR &&
                        fabs(x0.x() - tray->yhitList[iy]->pos.x()) < d2_width)  ynum++;
                        }
                }
        }
    return xnum+ynum;
}

int SiClusters::clusterSize (enum SiData::Axis a, unsigned int index)
{
  vector <Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  unsigned int nx = hitList->size();
  int size = 0;
  if(ix<nx) size = hitList->operator[](ix)->number;
  return size;
}

int SiClusters::tower(int indexhit)
{
        int itower = indexhit/1000000;
//      if (itower < 10)
//              std::cout << " warning : tower? "<< itower << "\n";
        return itower;
}

Point SiClusters::clusterPosition (enum SiData::Axis a, unsigned int index)
{
  vector <Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  unsigned int nx = hitList->size();
  //  int size = 0;
  Point Pos(0.,0.,0.);
  if(ix<nx) Pos = hitList->operator[](ix)->pos;
  return Pos;
}

int SiClusters::clusterNoise (enum SiData::Axis a, unsigned int index)
{
  vector <Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  if ( !(lyr < trayList.size() ) ) {
    // shouldn't happend, but it did
    return 0;
  }
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  if ( !hitList ) {
    // should never happen, but it did
    return 0;
  }
  unsigned int nx = hitList->size();
  int noise = 0;
  if(ix<nx) noise = hitList->operator[](ix)->num_noise;
  return noise;
}

bool SiClusters::hitFlagged (enum SiData::Axis a, unsigned int index)
{
  vector <Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  unsigned int nx = hitList->size();
  if(ix<nx) {
      if(hitList->operator[](ix)->flag) return true;
      else         return false;
  }
  return false;
}

void SiClusters::flagHit (unsigned int index)
{
  unsigned idx  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        break;
     }
  }
  if(tray) tray->flagHit(idx);
}

void SiClusters::flagHit (enum SiData::Axis a, unsigned int index)
{
  vector<Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  unsigned int nx = hitList->size();
  if(ix<nx) hitList->operator[](ix)->flag = 1;
}

void SiClusters::unflagHit (unsigned int index)
{
  unsigned idx  = index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        break;
     }
  }

  if(tray) tray->unflagHit(idx);
}

void SiClusters::unflagHit (enum SiData::Axis a, unsigned int index)
{
  vector<Cluster *> *hitList = 0;
  unsigned ix  =  index%10000;
  unsigned lyr  = (index/10000)%100;
  unsigned mod  = (index/1000000);

  Tray *tray = 0;
  vector<Tray *> &moduleTrays = trayList[lyr];
  for(unsigned j=0; j<moduleTrays.size(); j++) {
     if(moduleTrays[j]->module == mod) {
        tray = moduleTrays[j];
        if(a==SiData::X) hitList = &tray->xhitList;
        else             hitList = &tray->yhitList;

        break;
     }
  }
  unsigned int nx = hitList->size();
  if(ix<nx) hitList->operator[](ix)->flag = 0;
}

void SiClusters::clearHitFlags ()
{
  for(unsigned trayNum=0; trayNum < trayList.size(); trayNum++) {
     int nx= nHits(SiData::X,trayNum), ny= nHits(SiData::Y,trayNum);
     if( nx==0 && ny==0 ) continue;

     vector<Tray *> &moduleTrays = trayList[trayNum];
     unsigned j;
     for(j=0; j<moduleTrays.size(); j++)
         moduleTrays[j]->clearHitFlags();
   }
}

void SiClusters::printOn (std::ostream& cout)
{
  int trayNum, i, j;
  cout << "\nSiClusters:\n";
  for(trayNum=0; trayNum < (int)trayList.size(); trayNum++)
  {
     int nx= nHits(SiData::X,trayNum), ny= nHits(SiData::Y,trayNum);
     if( nx==0 && ny==0 ) continue;

     cout <<  trayNum << "X hits =" << nx <<" Y hits = " <<ny;
     vector<Tray *> &moduleTrays = trayList[trayNum];
     for(j=0; j< (int)moduleTrays.size(); j++) {
         Tray *tray = moduleTrays[j];
         unsigned modId = tray->module;
         nx = (int)tray->xhitList.size();
         cout <<  " Module " << modId << " X hits = " <<nx;
         for( i =0; i<nx; i++)
             cout << '\t'<< tray->xhitList[i]->pos << '\n';
         ny = (int)tray->yhitList.size();
         cout <<  " Module " << modId << " Y hits = " <<ny;
         for( i =0; i<ny; i++)
             cout << '\t'<< tray->yhitList[i]->pos << '\n';
      }
  }
  return;
}

SiClusters::Tray* SiClusters::getTrayFor (int lyr, int mod)
{
   // Does a Tray already exist for this layer/modules combination?
     vector<Tray *> &moduleTrays = trayList[lyr];
     for ( unsigned j=0; j<moduleTrays.size(); j++) {
       if ( (int)moduleTrays[j]->module == mod) {
           return moduleTrays[j];
       }
     }
    // Tray doesn't exist yet so make one...
      Tray *tray = new Tray(mod);
      trayList[lyr].push_back(tray);
      return tray;
}

void SiClusters::transferData (enum SiData::Axis a, const SiData& data,
                               int lyr, int start, int stop)
{
  // Form clusters of adjacent strips and store...
    float sum = 0., x, y, z;
    int number, noise =0;
    Point xSum(0., 0., 0.);
    unsigned lstIndex = 999999999;

    Tray *tray = getTrayFor(lyr, data.moduleId(a, lyr, start));

    for(int i=start; i<stop; i++){
        if(data.hitId(a, lyr, i) != lstIndex-1 &&
           data.hitId(a, lyr, i) != lstIndex+1 && lstIndex < 999999999) {

           number = (int)sum;
           x = xSum.x()/sum;
           y = xSum.y()/sum;
           z = xSum.z()/sum;
           if(a == SiData::X)

              tray->xhitList.push_back( new Cluster(Point(x,y,z), number, noise) );
           else
              tray->yhitList.push_back( new Cluster(Point(x,y,z), number, noise) );
           sum = 0.; noise = 0;
           xSum = Point(0., 0., 0.);
        }
        xSum += data.hit(a, lyr, i);
        sum  += 1.;
//      noise += data.hitType(a, lyr, i);
        lstIndex = data.hitId(a, lyr, i);
    }
    if(lstIndex< 999999999) {
       number = (int)sum;
       x = xSum.x()/sum;
       y = xSum.y()/sum;
       z = xSum.z()/sum;
       if(a == SiData::X)
          tray->xhitList.push_back( new Cluster(Point(x,y,z), number, noise) );
       else

          tray->yhitList.push_back( new Cluster(Point(x,y,z), number, noise) );
    }
}

---