---

mMuiFastRoadFinder.cxx

Go to the documentation of this file.

#include <iostream>
#include <math.h>

#include "PHNode.h"
#include "PHCompositeNode.h"
#include "PHDataNode.h"
#include "PHIODataNode.h"
#include "PHNodeIterator.h"

#include "mMuiFastRoadFinder.h"
#include "MuiGeomClasses.hh"

#include <fstream>

//MUIOO interfaces
#include <TMutNode.h>

// some LVL2 stuff; from lvl2_distribution/primitives
#include "Lvl2_MUIOO_MuiTriggerVars.hh"

static bool falseSignal = false;
static bool init_done = false;

const short NVIRT_PANELS = 8;

// Default constructor and destructor to pacify CINT
mMuiFastRoadFinder::mMuiFastRoadFinder() :
  _roadmap(0),
  _hitmap(0),
  _timer( PHTimeServer::get()->insert_new( "mMuiFastRoadFinder" ) )
{
  name = "mMuiFastRoadFinder";
}

PHBoolean 
mMuiFastRoadFinder::init()
{
  // Initialize signal array to zero.
  // ================================
  for (short logical = 0; logical < TMuiChannelId::kTwoPacksPerPlaneMax; logical++) {
    for (short gap = 0; gap < MUIOO::MAX_ARM; gap++) {
      _signal[logical][gap] = false;
    }
  }
  
  // Location pointed to by all symset members that don't have an 
  // actual tube to point to.
  // ============================================================
  for (short symset = 0; symset <  TMuiChannelId::kTwoPacksPerPlaneMax; 
       symset++) {
    for (short i = 0; i < mMuiFastRoadFinderPar::MAX_SYMSET; i++) {
      _trigger_signal[symset][i] = &falseSignal;
      short logical = symset + _mod_par->get_symset_logical_offset(i);
      if (logical >= 0 && 
          logical <  TMuiChannelId::kTwoPacksPerPlaneMax) {
        short gap = _mod_par->get_symset_plane(i);
        _trigger_signal[symset][i] = &(_signal[logical][gap]);
      }         
    }
  }
  // read the mapping files 
  // hardcode in default filenames for now
  read_idx_XoZ_map("/phenix/workarea/silvermy/offline/packages/mui/wrk/lvl2cosines.adb");
  read_twopack_idx_map("/phenix/workarea/silvermy/offline/packages/mui/wrk/muioomap.adb");
  return 0;
}

mMuiFastRoadFinder::~mMuiFastRoadFinder(){}

PHBoolean
mMuiFastRoadFinder::event(PHCompositeNode *top_node) 
{
  _timer.get()->restart();
  try {
    // Reset IOC pointers
    //
    set_interface_ptrs(top_node);

    if (!init_done)
      {
        init();
        init_done = true;
      }
    // Try to reconstruct some tracks
    //
    find_tracks();

  } catch(std::exception& e) {
    MUIOO::TRACE(e.what());
    return False;
  }
  _timer.get()->stop();
  return 0;
}

void 
mMuiFastRoadFinder::set_interface_ptrs(PHCompositeNode* top_node)
{
  // module runtime parameters
  _mod_par = TMutNode<mMuiFastRoadFinderPar>::find_node(top_node,"mMuiFastRoadFinderPar");
  // dump parameter listing
  //_mod_par->print();

  //Find the TMuiRoadO IOC
  _roadmap = TMutNode<TMuiRoadMapO>::find_node(top_node,"TMuiRoadMapO");
   
  //Find the TMuiHitO IOC
  _hitmap = TMutNode<TMuiHitMapO>::find_node(top_node,"TMuiHitMapO");
    
 return;  
}

void 
mMuiFastRoadFinder::read_idx_XoZ_map(const char* filename)
{
 // Read in Mui cosine/slope info from a file
  ifstream fin(filename);
  if (!fin)
    {
      cerr << PHWHERE << " - can't open " << filename
           << "for input" << endl;
      return;
    }  
  int nlines = 0;

  int idx = 0;
  int arm = 0;
  int orient = 0;
  float XoZ = 0.0;

  // read from the file as long as we can
  while ( fin.good() )
    {
      fin >> idx >> arm >> orient >> XoZ;

      if (! fin.good() ) 
        {
          break;
        }
        _idx_XoZ[arm][orient][idx] = XoZ;
        nlines++;
    }
  cout << PHWHERE << " XoZ nlines " << nlines << " read " << endl;
  return;
}


void 
mMuiFastRoadFinder::read_twopack_idx_map(const char* filename)
{
 // Read in Mui twopack->logical mapping from a file
  ifstream fin(filename);
  if (!fin)
    {
      cerr << PHWHERE << " - can't open " << filename
           << "for input" << endl;
      return;
    }  
  int nlines = 0;

  int arm = 0;
  int plane = 0;
  int panel = 0;
  int orient = 0;
  int twopack = 0;
  int idx = 0;


  // read from the file as long as we can
  while ( fin.good() )
    {
      fin >> arm >> plane >> panel >> orient >> twopack >> idx;

      if (! fin.good() ) 
        {
          break;
        }
        _twopack_idx[arm][plane][panel][orient][twopack] = idx;
        nlines++;
    }
  cout << PHWHERE << " twopack-idx nlines " << nlines << " read " << endl;
  return;
}

void
mMuiFastRoadFinder::find_tracks() 
{
  for (UShort_t arm = 0; arm < MUIOO::MAX_ARM; arm++)
    {
      for (UShort_t orient = 0; orient < MUIOO::MAX_ORIENTATION; orient++)
        {
          // diff in symset calc between original and simplified is hidden in 
          // a call in calc_symset method.
          _symsets[arm][orient].clear();

          if( _mod_par->get_mode() == mMuiFastRoadFinderPar::SIMPLIFIED ) {
            hit_to_logical(arm, orient); // fill signal variables
            calc_symset(arm, orient);
          }
          else if ( _mod_par->get_mode() == mMuiFastRoadFinderPar::ORIGINAL ) {
            // loop over virtual panels
            for (UShort_t ip = 0; ip<NVIRT_PANELS; ip++) {
              hit_to_logical_per_panel(arm, orient, ip); // fill signal variables
              calc_symset(arm, orient);
            }
          }
          else { // unknown mode - shouldn't be possible
            MUIOO::TRACE("mMuiFastRoadFinder::find_tracks - unknown alg. mode selected");
          }
        }

      //Make 2D tracks of all 1D combinations/symsets for this arm
      unsigned int hid = 0;
      unsigned int vid = 0;
      int depth = 0;
      mui_symset hsymset, vsymset;

      for(hid=0; hid<_symsets[arm][kHORIZ].size(); hid++) 
        {
          for(vid=0; vid<_symsets[arm][kVERT].size(); vid++) 
            {
              hsymset = _symsets[arm][kHORIZ][hid];
              vsymset = _symsets[arm][kVERT][vid];
              //Only construct if depth difference between orientations <= depthmatch
              if( abs(hsymset.depth - vsymset.depth) <= _mod_par->get_depth_match() )
                {
                  //The depth of the combination is the larger of the 1D depths
                  if (hsymset.depth > vsymset.depth)
                    {
                      depth = hsymset.depth;
                    }
                  else
                    {
                      depth = vsymset.depth;
                    }

                  // Note: maxmisshits cut and maxmisstotal cut based on gapbits are
                  // not included.

                  //Calculate direction of 2D Road. These are slopes not cosines
                  float cx1 = vsymset.cos;
                  float cy1 = hsymset.cos;
                  float invmag = 1.0/sqrt( cx1*cx1 + cy1*cy1 + 1 );

                  float dx = cx1*invmag; 
                  float dy = cy1*invmag; 
                  float dz = invmag;     

                  float slope = sqrt((dx*dx+dy*dy)/(dz*dz));
                  //float pmin = _mod_par->get_min_momentum(arm, depth)/invmag;

                  //Add the track to the list, if it passes cuts
                  if( depth >= _mod_par->get_min_depth() &&
                      slope >= _mod_par->get_min_slope() )
                    {
                      TMutFitPar fit2d;
                      TMuiRoadMapO::iterator iRoad = _roadmap->insert_new(arm);
                      TMuiRoadMapO::pointer cRoad = iRoad.current();

                      // Things to be set for real road:
                      // associated clusters are not really kept for now
                      
                      // nhit, depth, freedom, chi_square, gapbit 
                      cRoad->get()->set_nhit(hsymset.nhit + vsymset.nhit);
                                          
                      cRoad->get()->set_depth(depth);

                      cRoad->get()->set_freedom(hsymset.dof + vsymset.dof);

                      cRoad->get()->set_gapbit(hsymset.gapbit | 
                                               (vsymset.gapbit << MUIOO::MAX_PLANE));
                                          
                      //Combine the FitPar of the 1D roads                    
                      fit2d.set_x(vsymset.x);
                      fit2d.set_y(hsymset.y);
                      fit2d.set_z(hsymset.z);
                      fit2d.set_dxdz(cx1);
                      fit2d.set_dydz(cy1);
                      fit2d.set_z_end(hsymset.z);
                      fit2d.set_z_begin(hsymset.z);
                      Float_t chisq = (hsymset.chisq * hsymset.dof +
                                       vsymset.chisq * vsymset.dof)/
                        (hsymset.dof + vsymset.dof); 
                      fit2d.set_chi_square(chisq);

                      cRoad->get()->set_fit_par(fit2d);
                      // previously not set parameters
                      cRoad->get()->set_road_quality(chisq);
                      UShort_t gapmatch = hsymset.gapbit & 
                        vsymset.gapbit;
                      if (gapmatch != 0) {
                        cRoad->get()->set_max_hit_plane(2);
                      }
                      else {
                        cRoad->get()->set_max_hit_plane(1);
                      }
                    }
                } // depth_match
            } // vid
        } // hid
    } // arm

  return;
}

void 
mMuiFastRoadFinder::hit_to_logical(UShort_t arm, UShort_t orient)
{ // first clear the signal array
  memset(_signal,false,sizeof(_signal));
  
  UShort_t idx;
  
  TMuiHitMapO::iterator hit_iter = _hitmap->get(arm);
  while(TMuiHitMapO::pointer hit_ptr = hit_iter.next()) {
    if (hit_ptr->get()->get_orientation() == orient) {
      // right orientation selected
      
      // here we should get the 'logical tube'/twopack value
      // Check this again !!!
      idx = _twopack_idx[arm][hit_ptr->get()->get_plane()][hit_ptr->get()->get_panel()][orient][hit_ptr->get()->get_twopack()]; 

      _signal[idx][hit_ptr->get()->get_plane()] = true;
    }
  }
}

void 
mMuiFastRoadFinder::calc_symset(UShort_t arm, UShort_t orient)
{
  // Symset variables.
  // =================
  bool gap0bit[TMuiChannelId::kTwoPacksPerPlaneMax];
  bool gap1bit[TMuiChannelId::kTwoPacksPerPlaneMax];
  bool gap2bit[TMuiChannelId::kTwoPacksPerPlaneMax];
  bool gap3bit[TMuiChannelId::kTwoPacksPerPlaneMax];
  bool gap4bit[TMuiChannelId::kTwoPacksPerPlaneMax];
  
  // Loop over symsets; see if they contain a trigger road.
  // =====================================================
  for (short symset = 0; symset < TMuiChannelId::kTwoPacksPerPlaneMax; symset++) 
    {
      // Call the roadfinder, return bits that define tracking
      // status to each gap for each symset.
      // =====================================================
      gap0bit[symset]=false;
      gap1bit[symset]=false;
      gap2bit[symset]=false;
      gap3bit[symset]=false;
      gap4bit[symset]=false;
      
      //Only execute the road finder if there is a hit
      // in the (Note: middle of!) first two gaps
      if((*(_trigger_signal[symset][2]))
         ||(*(_trigger_signal[symset][10]))){
        bool** p = _trigger_signal[symset];
        gap0bit[symset]= *(_trigger_signal[symset][2]);

        if( _mod_par->get_mode() == mMuiFastRoadFinderPar::SIMPLIFIED ) {
          simple_trigger_road_finder(p, 
                                     gap1bit[symset], gap2bit[symset], 
                                     gap3bit[symset], gap4bit[symset]);
        }
        else if ( _mod_par->get_mode() == mMuiFastRoadFinderPar::ORIGINAL ) {
          orig_trigger_road_finder(p,
                                   gap1bit[symset], gap2bit[symset], 
                                   gap3bit[symset], gap4bit[symset]);
        }
        else { // unknown mode - shouldn't be possible
          MUIOO::TRACE("mMuiFastRoadFinder::calc_symset - unknown alg. mode selected");
        }
      } // one of 1st two gaps hit.
    } // symset 

  float idxXoZ[TMuiChannelId::kTwoPacksPerPlaneMax];
  for (short symset = 0; symset < TMuiChannelId::kTwoPacksPerPlaneMax; symset++) 
    {
      idxXoZ[symset] = _idx_XoZ[arm][orient][symset];
    }

  // Now let's summarize the results of the symset calculations
  // This includes rejecting adjacent, duplicate tracks
  mui_trigger_edge_finder(TMuiChannelId::kTwoPacksPerPlaneMax, idxXoZ,
                          gap0bit, gap1bit, gap2bit, gap3bit, gap4bit, 
                          arm, orient);
}

// Per-Panel stuff/Orig. alg.
void 
mMuiFastRoadFinder::hit_to_logical_per_panel(UShort_t arm, UShort_t orient,
                                             UShort_t virt_panel)
{ // first clear the signal array
  memset(_signal,false,sizeof(_signal));
  
  UShort_t idx;
  //We will populate hits per virtual panel which means. Panels:
  // 0 => Panel 0
  // 1 => Panel 0,1
  // 2 => Panel 1,2
  // 3 => Panel 2
  // 4 => Panel 3
  // 5 => Panel 3,4
  // 6 => Panel 4,5
  // 7 => Panel 5
  // In addition, we only want vertical tracks in panels 2,4 that have 
  //  slopes of appopriate sign.
  const short panel1[NVIRT_PANELS] = { 0, 1, 1, 2, 3, 4, 4, 5 };
  const short panel2[NVIRT_PANELS] = { 0, 0, 2, 2, 3, 3, 5, 5 };
  const float quadsign[NVIRT_PANELS] = { 1.0, 1.0, -1.0, -1.0, 
                                         -1.0, -1.0, 1.0, 1.0 };

  TMuiHitMapO::iterator hit_iter = _hitmap->get(arm);
  while(TMuiHitMapO::pointer hit_ptr = hit_iter.next()) {
    if (hit_ptr->get()->get_orientation() == orient) {
      // right orientation selected
      
      UShort_t hit_panel = hit_ptr->get()->get_panel();
      idx = _twopack_idx[arm][hit_ptr->get()->get_plane()][hit_ptr->get()->get_panel()][orient][hit_ptr->get()->get_twopack()]; 
      
      float hitslope = _idx_XoZ[arm][orient][idx];
      
      if( (hit_panel==panel1[virt_panel] || hit_panel==panel2[virt_panel]) 
          //Hit in virtual panel
         && ( (hit_ptr->get()->get_plane()>1) 
              // Accept all hits deeper than gap 1
              || (orient==1 && hitslope*quadsign[virt_panel]>0.0 &&
                  hit_panel==panel1[virt_panel])
              // Accept only Vertical with correct slope
              || (orient==0 && hit_panel==panel1[virt_panel])
              ) ) {
        //cout<<"Accepted\n";
        _signal[idx][hit_ptr->get()->get_plane()] = true;
      }
    }
  }
}

void 
mMuiFastRoadFinder::mui_trigger_edge_finder(int nsymsets, 
                    const float idx_XoZ[TMuiChannelId::kTwoPacksPerPlaneMax],
                                            bool *gap0bit,
                                            bool *gap1bit, bool *gap2bit, 
                                            bool *gap3bit, bool *gap4bit,
                                            UShort_t arm, UShort_t orient)
{
  // Zeroth elements have no previous bit to
  // negate them.
  // =======================================
  bool prev4 = gap4bit[0];
  bool prev3 = gap3bit[0];
  bool prev2 = gap2bit[0];
  bool prev1 = gap1bit[0];
  int i = 0;
  
  for (i = 1; i < nsymsets; i++) {
    bool curr4 = gap4bit[i];
    bool curr3 = gap3bit[i];
    bool curr2 = gap2bit[i];
    bool curr1 = gap1bit[i];
    gap4bit[i] = gap4bit[i] && !prev4;
    gap3bit[i] = gap3bit[i] && !prev3;
    gap2bit[i] = gap2bit[i] && !prev2;
    gap1bit[i] = gap1bit[i] && !prev1;
    prev4 = curr4;
    prev3 = curr3;
    prev2 = curr2;
    prev1 = curr1;
  }
  
  int nhit = 0;
  int gapbit = 0;

  // dummy values for now, for the degree of freedom and chisq
  int dof = 1;
  float chisq = 1.0; 
  // and for the coordinates
  float x = 0.0, y= 0.0, z= 0.0;

  for (i = 0; i < nsymsets; i++) {
    gapbit = 0;
    nhit = 0;
    if (gap0bit[i]) { gapbit |= 1; nhit++; }
    if (gap1bit[i]) { gapbit |= 2; nhit++; }
    if (gap2bit[i]) { gapbit |= 4; nhit++; }
    if (gap3bit[i]) { gapbit |= 8; nhit++; }
    if (gap4bit[i]) { gapbit |= 0x10; nhit++; }

    gap3bit[i] = gap3bit[i] || gap4bit[i];
    gap2bit[i] = gap2bit[i] || gap3bit[i];
    gap1bit[i] = gap1bit[i] || gap2bit[i];
    if(gap1bit[i])
      {
        mui_symset tmpSymset;
        tmpSymset.idx=i;
        tmpSymset.cos = idx_XoZ[i];
        tmpSymset.depth = 1;
        if(gap2bit[i])
          {
            tmpSymset.depth = 2;
          }
        if(gap3bit[i])
          {
            tmpSymset.depth = 3;
          }
        if(gap4bit[i])
          {
            tmpSymset.depth = 4;
          }
        // set all other values too
        tmpSymset.nhit= nhit;  
        tmpSymset.gapbit = gapbit;
        // even the dummies
        tmpSymset.dof = dof;   
        tmpSymset.chisq = chisq; 
        tmpSymset.x = x;     
        tmpSymset.y = y;     
        tmpSymset.z = z;     

        _symsets[arm][orient].push_back(tmpSymset);                
      }
  }
}

void 
mMuiFastRoadFinder::simple_trigger_road_finder(bool **p,
                                               bool& gap1bit, bool& gap2bit, 
                                               bool& gap3bit, bool& gap4bit)
{
  static const short search_width = 3;

  // these are counters for the indices in the gaps
  short i0 = 0, i1 = 0, i2 = 0, i3 = 0, i4 = 0;
  // gap 0:      0- 4; m = 2, 5 numbers
  // gap 1: 5  + 0-10; m = 10-5 = 5, 11 numbers
  // gap 2: 16 + 0-16; m = 24-16 = 8, 17 numbers
  // gap 3: 33 + 0-22; m = 44-33 = 11, 23 numbers
  // gap 4: 56 + 0-28; m = 70-56 = 14, 29 numbers
  // this info is stored in the mMuiFastRoadFinderPar.h file
  
  // negate the gapbits to start with
  gap1bit = false;
  gap2bit = false;
  gap3bit = false;
  gap4bit = false;
  short ngapsbest = 0;
  bool locgap0bit = false;
  bool locgap1bit = false;
  bool locgap2bit = false;
  bool locgap3bit = false;
  bool locgap4bit = false;
  bool candidateFound = false;
  
  //---------------------------------------------------------------------
  // The (software) algorithm works as follows:
  // loop over gap0 indices and see if we have something there or in gap1
  // for that index.
  // If so, we'll continue by looking at gaps 3 and 4.
  // we'll look X=3 wide in those last two gaps, i.e. -1,0,+1 + gap0-index
  // - i.e. use _mod_par->get_search_width() = 3; as default, to limit this
  // If we find something there also, we'll do the full OR between
  // -1,0,+1 + gap0-index in gaps 1 and 2 also. If the sum of these OR's
  // is larger or equal to 3, this is a validCandidate
  // - use const _mod_par->get_min_valid_hits() = 3;
  // we will then exit the search (gapbits should be filled properly)
  //
  // if we didn't find a valid candidate: set all gapbits to false
  //
  // New (Oct 31, 2003) implemented shallow check if no deep was found
  //---------------------------------------------------------------------
  
  i0 = 0; // start from the beginning
  short distFromMiddle; // help variable
  // start one below middle, by default
  short low_offset = 1;
  // go to one above middle
  // short high_offset = 1;
  // and search all
  short nsearch = search_width;
  int i;

  while ( i0<=_mod_par->get_symset_upper(0) )
    {
      distFromMiddle = i0 - _mod_par->get_symset_middle(0); // [-2,2]; 0 in the middle
      i1 = _mod_par->get_symset_middle(1) + distFromMiddle; // same index in gap1

      if ( *(p[i0]) || *(p[i1]) )
        { // ok, we have something in either gap 0 or 1
          // let's now get the OR's for gap 3 and 4.

          // Now, go through all the gaps..

          // gap 0
          locgap0bit = *(p[i0]);

          // gap 1
          i1 = _mod_par->get_symset_middle(1) + distFromMiddle - low_offset; 
          locgap1bit = false; 
          for (i = 0; i<nsearch; i++)
            {
              locgap1bit = locgap1bit || *(p[i1+i]);
            }

          // gap 2
          i2 = _mod_par->get_symset_middle(2) + distFromMiddle - low_offset; 
          locgap2bit = false; 
          for (i = 0; i<nsearch; i++)
            {
              locgap2bit = locgap2bit || *(p[i2+i]);
            }

          // gap 3
          i3 = _mod_par->get_symset_middle(3) + distFromMiddle - low_offset; 
          locgap3bit = false; 
          for (i = 0; i<nsearch; i++)
            {
              locgap3bit = locgap3bit || *(p[i3+i]);
            }

          // gap 4
          i4 = _mod_par->get_symset_middle(4) + distFromMiddle - low_offset; 
          locgap4bit = false;
          for (i = 0; i<nsearch; i++)
            {
              locgap4bit = locgap4bit || *(p[i4+i]);
            }

          // Requirement for sheep/deep:
          // if we have something either gap 3 or gap 4, we'll do the full OR
          // for gaps 1 and 2 also
          if (locgap3bit || locgap4bit)
            {
              short ngapsum = 0;
              if (locgap0bit) ngapsum++; 
              if (locgap1bit) ngapsum++;
              if (locgap2bit) ngapsum++;
              if (locgap3bit) ngapsum++;
              if (locgap4bit) ngapsum++;

              if (ngapsum >=_mod_par->get_min_valid_hits())
                {
                  candidateFound = true;
                  if (ngapsum > ngapsbest)
                    {
                      gap1bit = locgap1bit;
                      gap2bit = locgap2bit;
                      gap3bit = locgap3bit;
                      gap4bit = locgap4bit;
                    } // best so far
                } // ok
            } // fulfilled deep requirement
          else 
            { // could still perhaps be accepted as shallow
              // nothing in gap3 or 4, but we could perhaps still have enough
              // in gap0,1,2
              short ngapsum = 0;
              if (locgap0bit) ngapsum++; 
              if (locgap1bit) ngapsum++;
              if (locgap2bit) ngapsum++;
              if (ngapsum >=_mod_par->get_min_valid_hits())
                {
                  candidateFound = true;
                  if (ngapsum > ngapsbest)
                    {
                      gap1bit = locgap1bit;
                      gap2bit = locgap2bit;
                      gap3bit = locgap3bit;
                      gap4bit = locgap4bit;
                    } // best so far
                } // ok
            } // shallow test
        } // something in gap 0 or 1
      i0++;
    }

  // if we didn't find any candidate, we'll make sure that the
  // gapbits are all negated
  if (!candidateFound)
    {
      gap1bit = false;
      gap2bit = false;
      gap3bit = false;
      gap4bit = false;
    }

  // Return 4-bit value, the i^th bit indicates that a candidate made it
  // to the i^th gap
  // What we need to do is thus just go thru the lower gaps and see if
  // higher gaps were already set
  //
  // Outcommented this overwriting for now: DS, July 31, 2003 
  //
  /*
  gap3bit = gap3bit || gap4bit;
  gap2bit = gap2bit || gap3bit;
  gap1bit = gap1bit || gap2bit;
  */

  return;
}


void 
mMuiFastRoadFinder::orig_trigger_road_finder(bool **p,
                                             bool& gap1bit, bool& gap2bit, 
                                             bool& gap3bit, bool& gap4bit)
{
  short showermax[5][5]; // 1st index is maxdepth, 2nd is gap.
  bool roadfound[5][5];  // 1st index is maxdepth, 2nd is gap.
  // Search gap 0 for a shower (same code no matter the max depth).
  // ==============================================================
  gap0(p, &showermax[0][0]);
  roadfound[1][0] = *(p[m0]);
  roadfound[2][0] = *(p[m0]);
  roadfound[3][0] = *(p[m0]);
  roadfound[4][0] = *(p[m0]);
  // Search gap 1 for road, shower with maxdepth=1 cuts.
  // ===================================================
  gap1(p, _mod_par->get_swi(1,1), &showermax[1][1], &roadfound[1][1]);
          
  // Search gaps 1,2 for road, shower with maxdepth=2 cuts.
  // Gap 1 called again because of possible cut differences.
  // =======================================================
  gap1(p, _mod_par->get_swi(2,1), &showermax[2][1], &roadfound[2][1]);
  gap2(p, _mod_par->get_swi(2,2), &showermax[2][2], &roadfound[2][2]);
  
  // Search gaps 1,2,3 for road, shower with maxdepth=3 cuts.
  // Gap2 1,2 called again because of possible cut differences.
  // ==========================================================
  gap1(p, _mod_par->get_swi(3,1), &showermax[3][1], &roadfound[3][1]);
  gap2(p, _mod_par->get_swi(3,2), &showermax[3][2], &roadfound[3][2]);
  gap3(p, _mod_par->get_swi(3,3), &showermax[3][3], &roadfound[3][3]);
  
  // Search gaps 1-4 for road, shower with maxdepth=4 cuts.
  // Gaps 1-3 called again because of possible cut differences.
  // ==========================================================
  gap1(p, _mod_par->get_swi(4,1), &showermax[4][1], &roadfound[4][1]);
  gap2(p, _mod_par->get_swi(4,2), &showermax[4][2], &roadfound[4][2]);
  gap3(p, _mod_par->get_swi(4,3), &showermax[4][3], &roadfound[4][3]);
  gap4(p, _mod_par->get_swi(4,4), &showermax[4][4], &roadfound[4][4]);

  // We need to count the skipped gaps that were needed to generate 
  // each track candidate.
  // ==============================================================
  int nskip[5][5];      // 1st index is maxdepth, 2nd is gap.

  nskip[4][4] = 5 - bsum(&roadfound[4][0], 5);

  nskip[3][3] = 4 - bsum(&roadfound[3][0], 4);
  nskip[4][3] = 4 - bsum(&roadfound[4][0], 4);

  nskip[2][2] = 3 - bsum(&roadfound[2][0], 3);
  nskip[3][2] = 3 - bsum(&roadfound[3][0], 3);
  nskip[4][2] = 3 - bsum(&roadfound[4][0], 3);

  nskip[1][1] = 2 - bsum(&roadfound[1][0], 2);
  nskip[2][1] = 2 - bsum(&roadfound[2][0], 2);
  nskip[3][1] = 2 - bsum(&roadfound[3][0], 2);
  nskip[4][1] = 2 - bsum(&roadfound[4][0], 2);

  // Each maximum-depth=i candidate can find a road in any
  // gap j, for i>=j. Note that skipped gaps are handled in the 
  // assignment of *gapibits below.
  // ==========================================================
  bool road[5];

  road[4] = roadfound[4][4] && (nskip[4][4] <= _mod_par->get_maxskip(4,4));

  road[3] = (roadfound[3][3] && (nskip[3][3] <= _mod_par->get_maxskip(3,3))) ||
            (roadfound[4][3] && (nskip[4][3] <= _mod_par->get_maxskip(4,3))) ;
                                                         
  road[2] = (roadfound[2][2] && (nskip[2][2] <= _mod_par->get_maxskip(2,2))) ||
            (roadfound[3][2] && (nskip[3][2] <= _mod_par->get_maxskip(3,2))) ||
            (roadfound[4][2] && (nskip[4][2] <= _mod_par->get_maxskip(4,2))) ;
                                                         
  road[1] = (roadfound[1][1] && (nskip[1][1] <= _mod_par->get_maxskip(1,1))) ||
            (roadfound[2][1] && (nskip[2][1] <= _mod_par->get_maxskip(2,1))) ||
            (roadfound[3][1] && (nskip[3][1] <= _mod_par->get_maxskip(3,1))) ||
            (roadfound[4][1] && (nskip[4][1] <= _mod_par->get_maxskip(4,1)));

  // Conditions for setting showered bits:
  // 1. Kill any road candidate for which there is a more deeply
  //      reconstructed road that showers.
  // 2. Don't negate a more deeply reconstructed track, except by
  //      its own showermax cuts.
  // =============================================================
  bool showered[5];

  showered[4] = showermax[0][0] > _mod_par->get_showercut(4,0) ||
                showermax[4][1] > _mod_par->get_showercut(4,1) ||
                showermax[4][2] > _mod_par->get_showercut(4,2) ||
                showermax[4][3] > _mod_par->get_showercut(4,3) ||
                showermax[4][4] > _mod_par->get_showercut(4,4);

  showered[3] = showered[4] ||
                showermax[0][0] > _mod_par->get_showercut(3,0) ||
                showermax[3][1] > _mod_par->get_showercut(3,1) ||
                showermax[3][2] > _mod_par->get_showercut(3,2) ||
                showermax[3][3] > _mod_par->get_showercut(3,3);

  showered[2] = showered[3] ||
                showermax[0][0] > _mod_par->get_showercut(2,0) ||
                showermax[2][1] > _mod_par->get_showercut(2,1) ||
                showermax[2][2] > _mod_par->get_showercut(2,2);

  showered[1] = showered[2] ||
                showermax[0][0] > _mod_par->get_showercut(1,0) ||
                showermax[1][1] > _mod_par->get_showercut(1,1);

  // Return 4-bit value, the i^th bit indicates that a candidate made it
  // to the i^th gap w/o showering.
  // ===================================================================
  gap4bit =  road[4] && !showered[4];
  gap3bit = (road[3] && !showered[3]) || gap4bit;
  gap2bit = (road[2] && !showered[2]) || gap3bit;
  gap1bit = (road[1] && !showered[1]) || gap2bit;

  return;
}


void 
mMuiFastRoadFinder::gap0(bool **p, short *showermax)
{
  bool st[5];
  st[0] = *(p[k0]);
  st[1] = *(p[l0]);
  st[2] = *(p[m0]);
  st[3] = *(p[n0]);
  st[4] = *(p[o0]);
  
  *showermax = bsum(&st[0],5);
  m0c = *(p[m0]);
  m_0 = m0c;
 
  return;
}


void 
mMuiFastRoadFinder::gap1(bool **p, 
                         short swi, short *showermax, bool *roadfound)
{

  // These variables indicate whether a tube in the current gap is
  // pointed to by a tube in the previous gap. 
  // Note: The way this is done is repetitious for gap1, but
  //       I did it because it made it easier
  //       for me to get the pattern for subsequent gaps.
  // =============================================================
  bool j1m0 = m_0 && swi==3; bool j1m0c = m0c && swi==3;
  bool k1m0 = m_0 && swi>=2; bool k1m0c = m0c && swi>=2;
  bool l1m0 = m_0;           bool l1m0c = m0c;            
  bool m1m0 = m_0;           bool m1m0c = m0c;            
  bool n1m0 = m_0;           bool n1m0c = m0c;            
  bool o1m0 = m_0 && swi>=2; bool o1m0c = m0c && swi>=2;
  bool p1m0 = m_0 && swi==3; bool p1m0c = m0c && swi==3;

  bool j1pre = j1m0; 
  bool k1pre = k1m0; 
  bool l1pre = l1m0;   
  bool m1pre = true ; // Allow skipped gap0.   
  bool n1pre = n1m0;   
  bool o1pre = o1m0; 
  bool p1pre = p1m0; 
  
  // Search for road in current gap.
  // ===============================
  j1c = *(p[musj1]) && j1pre;
  k1c = *(p[k1]) && k1pre;
  l1c = *(p[l1]) && l1pre;
  m1c = *(p[m1]) && m1pre;
  n1c = *(p[n1]) && n1pre;
  o1c = *(p[o1]) && o1pre;
  p1c = *(p[p1]) && p1pre;
  
  *roadfound = j1c || k1c || l1c || m1c 
                   || n1c || o1c || p1c;
  
  // Find shower size.
  // =================
  static const short nInShower = 5;
  static const short nShowerSums = 7;
  bool st[nShowerSums+nInShower-1];
  st[0]  = *(p[h1]) && j1pre;
  st[1]  = *(p[i1]) && k1pre;
  st[2]  = *(p[musj1]) && l1pre;
  st[3]  = *(p[k1]) && m1pre;
  st[4]  = *(p[l1]) && m1pre;
  st[5]  = *(p[m1]) && m1pre;
  st[6]  = *(p[n1]) && m1pre;
  st[7]  = *(p[o1]) && m1pre;
  st[8]  = *(p[p1]) && n1pre;
  st[9]  = *(p[q1]) && o1pre;
  st[10] = *(p[r1]) && p1pre;
  
  *showermax = 0;
  for (short i = 0; i < nShowerSums; i++) {
    short shower = (short) bsum(&st[i],nInShower);
    if (shower > *showermax)
      *showermax = shower;
  }

  // We need to find out if this gap was skipped.
  // The first step is to define variables which indicate
  // whether the set tubes on this gap which are projected 
  // to by each tube on the previous gap are MT.
  // Note: The way this is done is repetitious for gap1, but
  //       I did it because it made it easier
  //       for me to get the pattern for subsequent gaps.
  // ==========================================================
  bool m0sl = *(p[l1]) || *(p[m1]) || *(p[n1]) ||
    (swi == 2 && (*(p[k1]) || *(p[o1]))) ||
    (swi == 3 && (*(p[musj1]) || *(p[p1])));

  // If a tube on the current gap is pointed to by a tube 
  // the previous gap, then the current-gap tube's signal projects forward to 
  // the next gap iff all tubes pointed to by the previous-gap tube are MT.
  // ========================================================================
  bool j1m0skip1 = j1m0c && !m0sl;
  bool k1m0skip1 = k1m0c && !m0sl;
  bool l1m0skip1 = l1m0c && !m0sl;
  bool m1m0skip1 = m1m0c && !m0sl;
  bool n1m0skip1 = n1m0c && !m0sl;
  bool o1m0skip1 = o1m0c && !m0sl;
  bool p1m0skip1 = p1m0c && !m0sl;

  bool j1skip1 = j1m0skip1;
  bool k1skip1 = k1m0skip1;
  bool l1skip1 = l1m0skip1;
  bool m1skip1 = m1m0skip1;
  bool n1skip1 = n1m0skip1;
  bool o1skip1 = o1m0skip1;
  bool p1skip1 = p1m0skip1;

  // Projection to the next gap.
  // ===========================
  j_1 = j1c || j1skip1;
  k_1 = k1c || k1skip1;
  l_1 = l1c || l1skip1;
  m_1 = m1c || m1skip1;
  n_1 = n1c || n1skip1;
  o_1 = o1c || o1skip1;
  p_1 = p1c || p1skip1;
  
  return;
}


void 
mMuiFastRoadFinder::gap2(bool **p, 
                         short swi, short *showermax, bool *roadfound)
{

  // These variables indicate whether a tube in this gap is
  // pointed to by a tube in the previous gap. 
  // ======================================================
  bool g2j1 = j_1 && swi==3;   bool g2j1c = j1c && swi==3;
  bool h2j1 = j_1 && swi>=2;   bool h2j1c = j1c && swi>=2;
  bool i2j1 = j_1;             bool i2j1c = j1c;            
  bool j2j1 = j_1;             bool j2j1c = j1c;            
  bool k2j1 = j_1;             bool k2j1c = j1c;            
  bool l2j1 = j_1 && swi>=2;   bool l2j1c = j1c && swi>=2;
  bool m2j1 = j_1 && swi==3;   bool m2j1c = j1c && swi==3;
                                                         
  bool h2k1 = k_1 && swi==3;   bool h2k1c = k1c && swi==3;
  bool i2k1 = k_1 && swi>=2;   bool i2k1c = k1c && swi>=2;
  bool j2k1 = k_1;             bool j2k1c = k1c;            
  bool k2k1 = k_1;             bool k2k1c = k1c;            
  bool l2k1 = k_1;             bool l2k1c = k1c;            
  bool m2k1 = k_1 && swi>=2;   bool m2k1c = k1c && swi>=2;
  bool n2k1 = k_1 && swi==3;   bool n2k1c = k1c && swi==3;
                                                         
  bool i2l1 = l_1 && swi==3;   bool i2l1c = l1c && swi==3;
  bool j2l1 = l_1 && swi>=2;   bool j2l1c = l1c && swi>=2;
  bool k2l1 = l_1;             bool k2l1c = l1c;            
  bool l2l1 = l_1;             bool l2l1c = l1c;            
  bool m2l1 = l_1;             bool m2l1c = l1c;            
  bool n2l1 = l_1 && swi>=2;   bool n2l1c = l1c && swi>=2;
  bool o2l1 = l_1 && swi==3;   bool o2l1c = l1c && swi==3;
                                                         
  bool j2m1 = m_1 && swi==3;   bool j2m1c = m1c && swi==3;
  bool k2m1 = m_1 && swi>=2;   bool k2m1c = m1c && swi>=2;
  bool l2m1 = m_1;             bool l2m1c = m1c;            
  bool m2m1 = m_1;             bool m2m1c = m1c;            
  bool n2m1 = m_1;             bool n2m1c = m1c;            
  bool o2m1 = m_1 && swi>=2;   bool o2m1c = m1c && swi>=2;
  bool p2m1 = m_1 && swi==3;   bool p2m1c = m1c && swi==3;
                                                         
  bool k2n1 = n_1 && swi==3;   bool k2n1c = n1c && swi==3;
  bool l2n1 = n_1 && swi>=2;   bool l2n1c = n1c && swi>=2;
  bool m2n1 = n_1;             bool m2n1c = n1c;            
  bool n2n1 = n_1;             bool n2n1c = n1c;            
  bool o2n1 = n_1;             bool o2n1c = n1c;            
  bool p2n1 = n_1 && swi>=2;   bool p2n1c = n1c && swi>=2;
  bool q2n1 = n_1 && swi==3;   bool q2n1c = n1c && swi==3;
                                                         
  bool l2o1 = o_1 && swi==3;   bool l2o1c = o1c && swi==3;
  bool m2o1 = o_1 && swi>=2;   bool m2o1c = o1c && swi>=2;
  bool n2o1 = o_1;             bool n2o1c = o1c;            
  bool o2o1 = o_1;             bool o2o1c = o1c;            
  bool p2o1 = o_1;             bool p2o1c = o1c;            
  bool q2o1 = o_1 && swi>=2;   bool q2o1c = o1c && swi>=2;
  bool r2o1 = o_1 && swi==3;   bool r2o1c = o1c && swi==3;
                                                         
  bool m2p1 = p_1 && swi==3;   bool m2p1c = p1c && swi==3;
  bool n2p1 = p_1 && swi>=2;   bool n2p1c = p1c && swi>=2;
  bool o2p1 = p_1;             bool o2p1c = p1c;            
  bool p2p1 = p_1;             bool p2p1c = p1c;            
  bool q2p1 = p_1;             bool q2p1c = p1c;            
  bool r2p1 = p_1 && swi>=2;   bool r2p1c = p1c && swi>=2;
  bool s2p1 = p_1 && swi==3;   bool s2p1c = p1c && swi==3;

  bool g2pre = g2j1; 
  bool h2pre = h2j1 || h2k1;
  bool i2pre = i2j1 || i2k1 || i2l1;
  bool j2pre = j2j1 || j2k1 || j2l1 || j2m1;
  bool k2pre = k2j1 || k2k1 || k2l1 || k2m1 || k2n1;
  bool l2pre = l2j1 || l2k1 || l2l1 || l2m1 || l2n1 || l2o1;
  bool m2pre = m2j1 || m2k1 || m2l1 || m2m1 || m2n1 || m2o1 || m2p1;
  bool n2pre = n2k1 || n2l1 || n2m1 || n2n1 || n2o1 || n2p1;
  bool o2pre = o2l1 || o2m1 || o2n1 || o2o1 || o2p1;       
  bool p2pre = p2m1 || p2n1 || p2o1 || p2p1;
  bool q2pre = q2n1 || q2o1 || q2p1;
  bool r2pre = r2o1 || r2p1;
  bool s2pre = s2p1;
  
  // Search for road in current gap.
  // ===============================
  g2c = *(p[g2]) && g2pre;
  h2c = *(p[h2]) && h2pre;
  i2c = *(p[i2]) && i2pre;
  j2c = *(p[j2]) && j2pre;
  k2c = *(p[k2]) && k2pre;
  l2c = *(p[l2]) && l2pre;
  m2c = *(p[m2]) && m2pre;
  n2c = *(p[n2]) && n2pre;
  o2c = *(p[o2]) && o2pre;
  p2c = *(p[p2]) && p2pre;
  q2c = *(p[q2]) && q2pre;
  r2c = *(p[r2]) && r2pre;
  s2c = *(p[s2]) && s2pre;
  
  *roadfound = g2c || h2c || i2c || j2c || k2c || l2c || m2c 
                   || n2c || o2c || p2c || q2c || r2c || s2c;
  
  // Find shower size.
  // =================
  static const short nInShower = 5;
  static const short nShowerSums = 13;
  bool st[nShowerSums+nInShower-1];
  st[0]  = *(p[e2]) && g2pre;
  st[1]  = *(p[f2]) && h2pre;
  st[2]  = *(p[g2]) && i2pre;
  st[3]  = *(p[h2]) && j2pre;
  st[4]  = *(p[i2]) && k2pre;
  st[5]  = *(p[j2]) && l2pre;
  st[6]  = *(p[k2]) && m2pre;
  st[7]  = *(p[l2]) && m2pre;
  st[8]  = *(p[m2]) && m2pre;
  st[9]  = *(p[n2]) && m2pre;
  st[10] = *(p[o2]) && m2pre;
  st[11] = *(p[p2]) && n2pre;
  st[12] = *(p[q2]) && o2pre;
  st[13] = *(p[r2]) && p2pre;
  st[14] = *(p[s2]) && q2pre;
  st[15] = *(p[t2]) && r2pre;
  st[16] = *(p[u2]) && s2pre;
  
  *showermax = 0;
  for (short i = 0; i < nShowerSums; i++) {
    short shower = (short) bsum(&st[i],nInShower);
    if (shower > *showermax)
      *showermax = shower;
  }

  // We need to find out if this gap was skipped.
  // The first step is to define variables which indicate
  // whether the set tubes on this gap which are projected 
  // to by each tube on the previous gap are MT.
  // =====================================================
  bool j1sl = *(p[i2]) || *(p[j2]) || *(p[k2]) ||
    (swi == 2 && (*(p[h2]) || *(p[l2]))) ||
    (swi == 3 && (*(p[g2]) || *(p[m2])));

  bool k1sl = *(p[j2]) || *(p[k2]) || *(p[l2]) ||
    (swi == 2 && (*(p[i2]) || *(p[m2]))) ||
    (swi == 3 && (*(p[h2]) || *(p[n2])));

  bool l1sl = *(p[k2]) || *(p[l2]) || *(p[m2]) ||
    (swi == 2 && (*(p[j2]) || *(p[n2]))) ||
    (swi == 3 && (*(p[i2]) || *(p[o2])));

  bool m1sl = *(p[l2]) || *(p[m2]) || *(p[n2]) ||
    (swi == 2 && (*(p[k2]) || *(p[o2]))) ||
    (swi == 3 && (*(p[j2]) || *(p[p2])));

  bool n1sl = *(p[m2]) || *(p[n2]) || *(p[o2]) ||
    (swi == 2 && (*(p[l2]) || *(p[p2]))) ||
    (swi == 3 && (*(p[k2]) || *(p[q2])));

  bool o1sl = *(p[n2]) || *(p[o2]) || *(p[p2]) ||
    (swi == 2 && (*(p[m2]) || *(p[q2]))) ||
    (swi == 3 && (*(p[l2]) || *(p[r2])));

  bool p1sl = *(p[o2]) || *(p[p2]) || *(p[q2]) ||
    (swi == 2 && (*(p[n2]) || *(p[r2]))) ||
    (swi == 3 && (*(p[m2]) || *(p[s2])));

  // If a tube on the current gap is pointed to by a tube 
  // the previous gap, then the current-gap tube's signal projects forward to 
  // the next gap iff all tubes pointed to by the previous-gap tube are MT.
  // ========================================================================
  bool g2j1skip2 = g2j1c && !j1sl;
  bool h2j1skip2 = h2j1c && !j1sl;
  bool i2j1skip2 = i2j1c && !j1sl;
  bool j2j1skip2 = j2j1c && !j1sl;
  bool k2j1skip2 = k2j1c && !j1sl;
  bool l2j1skip2 = l2j1c && !j1sl;
  bool m2j1skip2 = m2j1c && !j1sl;
                              
  bool h2k1skip2 = h2k1c && !k1sl;
  bool i2k1skip2 = i2k1c && !k1sl;
  bool j2k1skip2 = j2k1c && !k1sl;
  bool k2k1skip2 = k2k1c && !k1sl;
  bool l2k1skip2 = l2k1c && !k1sl;
  bool m2k1skip2 = m2k1c && !k1sl;
  bool n2k1skip2 = n2k1c && !k1sl;
                              
  bool i2l1skip2 = i2l1c && !l1sl;
  bool j2l1skip2 = j2l1c && !l1sl;
  bool k2l1skip2 = k2l1c && !l1sl;
  bool l2l1skip2 = l2l1c && !l1sl;
  bool m2l1skip2 = m2l1c && !l1sl;
  bool n2l1skip2 = n2l1c && !l1sl;
  bool o2l1skip2 = o2l1c && !l1sl;
                              
  bool j2m1skip2 = j2m1c && !m1sl;
  bool k2m1skip2 = k2m1c && !m1sl;
  bool l2m1skip2 = l2m1c && !m1sl;
  bool m2m1skip2 = m2m1c && !m1sl;
  bool n2m1skip2 = n2m1c && !m1sl;
  bool o2m1skip2 = o2m1c && !m1sl;
  bool p2m1skip2 = p2m1c && !m1sl;
                              
  bool k2n1skip2 = k2n1c && !n1sl;
  bool l2n1skip2 = l2n1c && !n1sl;
  bool m2n1skip2 = m2n1c && !n1sl;
  bool n2n1skip2 = n2n1c && !n1sl;
  bool o2n1skip2 = o2n1c && !n1sl;
  bool p2n1skip2 = p2n1c && !n1sl;
  bool q2n1skip2 = q2n1c && !n1sl;
                              
  bool l2o1skip2 = l2o1c && !o1sl;
  bool m2o1skip2 = m2o1c && !o1sl;
  bool n2o1skip2 = n2o1c && !o1sl;
  bool o2o1skip2 = o2o1c && !o1sl;
  bool p2o1skip2 = p2o1c && !o1sl;
  bool q2o1skip2 = q2o1c && !o1sl;
  bool r2o1skip2 = r2o1c && !o1sl;
                              
  bool m2p1skip2 = m2p1c && !p1sl;
  bool n2p1skip2 = n2p1c && !p1sl;
  bool o2p1skip2 = o2p1c && !p1sl;
  bool p2p1skip2 = p2p1c && !p1sl;
  bool q2p1skip2 = q2p1c && !p1sl;
  bool r2p1skip2 = r2p1c && !p1sl;
  bool s2p1skip2 = s2p1c && !p1sl;

  bool g2skip2 = g2j1skip2;
  bool h2skip2 = h2j1skip2 || h2k1skip2;
  bool i2skip2 = i2j1skip2 || i2k1skip2 || i2l1skip2;
  bool j2skip2 = j2j1skip2 || j2k1skip2 || j2l1skip2 || j2m1skip2;
  bool k2skip2 = k2j1skip2 || k2k1skip2 || k2l1skip2 || k2m1skip2 || k2n1skip2;
  bool l2skip2 = l2j1skip2 || l2k1skip2 || l2l1skip2 || l2m1skip2 || l2n1skip2 || l2o1skip2;
  bool m2skip2 = m2j1skip2 || m2k1skip2 || m2l1skip2 || m2m1skip2 || m2n1skip2 || m2o1skip2 || m2p1skip2;
  bool n2skip2 = n2k1skip2 || n2l1skip2 || n2m1skip2 || n2n1skip2 || n2o1skip2 || n2p1skip2;
  bool o2skip2 = o2l1skip2 || o2m1skip2 || o2n1skip2 || o2o1skip2 || o2p1skip2;       
  bool p2skip2 = p2m1skip2 || p2n1skip2 || p2o1skip2 || p2p1skip2;
  bool q2skip2 = q2n1skip2 || q2o1skip2 || q2p1skip2;
  bool r2skip2 = r2o1skip2 || r2p1skip2;
  bool s2skip2 = s2p1skip2;
  
  // Projection to the next gap.
  // ===========================
  g_2 = g2c || g2skip2;
  h_2 = h2c || h2skip2;
  i_2 = i2c || i2skip2;
  j_2 = j2c || j2skip2;
  k_2 = k2c || k2skip2;
  l_2 = l2c || l2skip2;
  m_2 = m2c || m2skip2;
  n_2 = n2c || n2skip2;
  o_2 = o2c || o2skip2;
  p_2 = p2c || p2skip2;
  q_2 = q2c || q2skip2;
  r_2 = r2c || r2skip2;
  s_2 = s2c || s2skip2;
  
  return;
}



void 
mMuiFastRoadFinder::gap3(bool **p, 
                         short swi, short *showermax, bool *roadfound)
{

  // These variables indicate whether a tube in this gap is
  // pointed to by a tube in the previous gap. 
  // ======================================================
  bool d3g2 = g_2 && swi==3;   bool d3g2c = g2c && swi==3;
  bool e3g2 = g_2 && swi>=2;   bool e3g2c = g2c && swi>=2;
  bool f3g2 = g_2;             bool f3g2c = g2c;            
  bool g3g2 = g_2;             bool g3g2c = g2c;            
  bool h3g2 = g_2;             bool h3g2c = g2c;            
  bool i3g2 = g_2 && swi>=2;   bool i3g2c = g2c && swi>=2;
  bool j3g2 = g_2 && swi==3;   bool j3g2c = g2c && swi==3;
                                                         
  bool e3h2 = h_2 && swi==3;   bool e3h2c = h2c && swi==3;
  bool f3h2 = h_2 && swi>=2;   bool f3h2c = h2c && swi>=2;
  bool g3h2 = h_2;             bool g3h2c = h2c;            
  bool h3h2 = h_2;             bool h3h2c = h2c;            
  bool i3h2 = h_2;             bool i3h2c = h2c;            
  bool j3h2 = h_2 && swi>=2;   bool j3h2c = h2c && swi>=2;
  bool k3h2 = h_2 && swi==3;   bool k3h2c = h2c && swi==3;
                                                         
  bool f3i2 = i_2 && swi==3;   bool f3i2c = i2c && swi==3;
  bool g3i2 = i_2 && swi>=2;   bool g3i2c = i2c && swi>=2;
  bool h3i2 = i_2;             bool h3i2c = i2c;            
  bool i3i2 = i_2;             bool i3i2c = i2c;            
  bool j3i2 = i_2;             bool j3i2c = i2c;            
  bool k3i2 = i_2 && swi>=2;   bool k3i2c = i2c && swi>=2;
  bool l3i2 = i_2 && swi==3;   bool l3i2c = i2c && swi==3;
                                                         
  bool g3j2 = j_2 && swi==3;   bool g3j2c = j2c && swi==3;
  bool h3j2 = j_2 && swi>=2;   bool h3j2c = j2c && swi>=2;
  bool i3j2 = j_2;             bool i3j2c = j2c;            
  bool j3j2 = j_2;             bool j3j2c = j2c;            
  bool k3j2 = j_2;             bool k3j2c = j2c;            
  bool l3j2 = j_2 && swi>=2;   bool l3j2c = j2c && swi>=2;
  bool m3j2 = j_2 && swi==3;   bool m3j2c = j2c && swi==3;
                                                         
  bool h3k2 = k_2 && swi==3;   bool h3k2c = k2c && swi==3;
  bool i3k2 = k_2 && swi>=2;   bool i3k2c = k2c && swi>=2;
  bool j3k2 = k_2;             bool j3k2c = k2c;            
  bool k3k2 = k_2;             bool k3k2c = k2c;            
  bool l3k2 = k_2;             bool l3k2c = k2c;            
  bool m3k2 = k_2 && swi>=2;   bool m3k2c = k2c && swi>=2;
  bool n3k2 = k_2 && swi==3;   bool n3k2c = k2c && swi==3;
                                                         
  bool i3l2 = l_2 && swi==3;   bool i3l2c = l2c && swi==3;
  bool j3l2 = l_2 && swi>=2;   bool j3l2c = l2c && swi>=2;
  bool k3l2 = l_2;             bool k3l2c = l2c;            
  bool l3l2 = l_2;             bool l3l2c = l2c;            
  bool m3l2 = l_2;             bool m3l2c = l2c;            
  bool n3l2 = l_2 && swi>=2;   bool n3l2c = l2c && swi>=2;
  bool o3l2 = l_2 && swi==3;   bool o3l2c = l2c && swi==3;
                                                         
  bool j3m2 = m_2 && swi==3;   bool j3m2c = m2c && swi==3;
  bool k3m2 = m_2 && swi>=2;   bool k3m2c = m2c && swi>=2;
  bool l3m2 = m_2;             bool l3m2c = m2c;            
  bool m3m2 = m_2;             bool m3m2c = m2c;            
  bool n3m2 = m_2;             bool n3m2c = m2c;            
  bool o3m2 = m_2 && swi>=2;   bool o3m2c = m2c && swi>=2;
  bool p3m2 = m_2 && swi==3;   bool p3m2c = m2c && swi==3;
                                                         
  bool k3n2 = n_2 && swi==3;   bool k3n2c = n2c && swi==3;
  bool l3n2 = n_2 && swi>=2;   bool l3n2c = n2c && swi>=2;
  bool m3n2 = n_2;             bool m3n2c = n2c;            
  bool n3n2 = n_2;             bool n3n2c = n2c;            
  bool o3n2 = n_2;             bool o3n2c = n2c;            
  bool p3n2 = n_2 && swi>=2;   bool p3n2c = n2c && swi>=2;
  bool q3n2 = n_2 && swi==3;   bool q3n2c = n2c && swi==3;
                                                         
  bool l3o2 = o_2 && swi==3;   bool l3o2c = o2c && swi==3;
  bool m3o2 = o_2 && swi>=2;   bool m3o2c = o2c && swi>=2;
  bool n3o2 = o_2;             bool n3o2c = o2c;            
  bool o3o2 = o_2;             bool o3o2c = o2c;            
  bool p3o2 = o_2;             bool p3o2c = o2c;            
  bool q3o2 = o_2 && swi>=2;   bool q3o2c = o2c && swi>=2;
  bool r3o2 = o_2 && swi==3;   bool r3o2c = o2c && swi==3;
                                                         
  bool m3p2 = p_2 && swi==3;   bool m3p2c = p2c && swi==3;
  bool n3p2 = p_2 && swi>=2;   bool n3p2c = p2c && swi>=2;
  bool o3p2 = p_2;             bool o3p2c = p2c;            
  bool p3p2 = p_2;             bool p3p2c = p2c;            
  bool q3p2 = p_2;             bool q3p2c = p2c;            
  bool r3p2 = p_2 && swi>=2;   bool r3p2c = p2c && swi>=2;
  bool s3p2 = p_2 && swi==3;   bool s3p2c = p2c && swi==3;
                                                         
  bool n3q2 = q_2 && swi==3;   bool n3q2c = q2c && swi==3;
  bool o3q2 = q_2 && swi>=2;   bool o3q2c = q2c && swi>=2;
  bool p3q2 = q_2;             bool p3q2c = q2c;            
  bool q3q2 = q_2;             bool q3q2c = q2c;            
  bool r3q2 = q_2;             bool r3q2c = q2c;            
  bool s3q2 = q_2 && swi>=2;   bool s3q2c = q2c && swi>=2;
  bool t3q2 = q_2 && swi==3;   bool t3q2c = q2c && swi==3;
                                                         
  bool o3r2 = r_2 && swi==3;   bool o3r2c = r2c && swi==3;
  bool p3r2 = r_2 && swi>=2;   bool p3r2c = r2c && swi>=2;
  bool q3r2 = r_2;             bool q3r2c = r2c;            
  bool r3r2 = r_2;             bool r3r2c = r2c;            
  bool s3r2 = r_2;             bool s3r2c = r2c;            
  bool t3r2 = r_2 && swi>=2;   bool t3r2c = r2c && swi>=2;
  bool u3r2 = r_2 && swi==3;   bool u3r2c = r2c && swi==3;
                                                         
  bool p3s2 = s_2 && swi==3;   bool p3s2c = s2c && swi==3;
  bool q3s2 = s_2 && swi>=2;   bool q3s2c = s2c && swi>=2;
  bool r3s2 = s_2;             bool r3s2c = s2c;            
  bool s3s2 = s_2;             bool s3s2c = s2c;            
  bool t3s2 = s_2;             bool t3s2c = s2c;            
  bool u3s2 = s_2 && swi>=2;   bool u3s2c = s2c && swi>=2;
  bool v3s2 = s_2 && swi==3;   bool v3s2c = s2c && swi==3;

  bool d3pre = d3g2; 
  bool e3pre = e3g2 || e3h2;
  bool f3pre = f3g2 || f3h2 || f3i2;
  bool g3pre = g3g2 || g3h2 || g3i2 || g3j2;
  bool h3pre = h3g2 || h3h2 || h3i2 || h3j2 || h3k2;
  bool i3pre = i3g2 || i3h2 || i3i2 || i3j2 || i3k2 || i3l2;
  bool j3pre = j3g2 || j3h2 || j3i2 || j3j2 || j3k2 || j3l2 || j3m2;
  bool k3pre = k3h2 || k3i2 || k3j2 || k3k2 || k3l2 || k3m2 || k3n2;
  bool l3pre = l3i2 || l3j2 || l3k2 || l3l2 || l3m2 || l3n2 || l3o2;
  bool m3pre = m3j2 || m3k2 || m3l2 || m3m2 || m3n2 || m3o2 || m3p2;
  bool n3pre = n3k2 || n3l2 || n3m2 || n3n2 || n3o2 || n3p2 || n3q2;
  bool o3pre = o3l2 || o3m2 || o3n2 || o3o2 || o3p2 || o3q2 || o3r2;
  bool p3pre = p3m2 || p3n2 || p3o2 || p3p2 || p3q2 || p3r2 || p3s2;
  bool q3pre = q3n2 || q3o2 || q3p2 || q3q2 || q3r2 || q3s2;
  bool r3pre = r3o2 || r3p2 || r3q2 || r3r2 || r3s2;       
  bool s3pre = s3p2 || s3q2 || s3r2 || s3s2;
  bool t3pre = t3q2 || t3r2 || t3s2;
  bool u3pre = u3r2 || u3s2;
  bool v3pre = v3s2;
  
  // Search for road in current gap.
  // ===============================
  d3c = *(p[d3]) && d3pre;
  e3c = *(p[e3]) && e3pre;
  f3c = *(p[f3]) && f3pre;
  g3c = *(p[g3]) && g3pre;
  h3c = *(p[h3]) && h3pre;
  i3c = *(p[i3]) && i3pre;
  j3c = *(p[j3]) && j3pre;
  k3c = *(p[k3]) && k3pre;
  l3c = *(p[l3]) && l3pre;
  m3c = *(p[m3]) && m3pre;
  n3c = *(p[n3]) && n3pre;
  o3c = *(p[o3]) && o3pre;
  p3c = *(p[p3]) && p3pre;
  q3c = *(p[q3]) && q3pre;
  r3c = *(p[r3]) && r3pre;
  s3c = *(p[s3]) && s3pre;
  t3c = *(p[t3]) && t3pre;
  u3c = *(p[u3]) && u3pre;
  v3c = *(p[v3]) && v3pre;

  *roadfound = d3c || e3c || f3c || g3c || h3c || i3c || j3c || k3c || l3c || m3c 
                   || n3c || o3c || p3c || q3c || r3c || s3c || t3c || u3c || v3c;
  
  // Find shower size.
  // =================
  static const short nInShower = 5;
  static const short nShowerSums = 19;
  bool st[nShowerSums+nInShower-1];
  st[0]  = *(p[b3]) && d3pre;
  st[1]  = *(p[c3]) && e3pre;
  st[2]  = *(p[d3]) && f3pre;
  st[3]  = *(p[e3]) && g3pre;
  st[4]  = *(p[f3]) && h3pre;
  st[5]  = *(p[g3]) && i3pre;
  st[6]  = *(p[h3]) && j3pre;
  st[7]  = *(p[i3]) && k3pre;
  st[8]  = *(p[j3]) && l3pre;
  st[9]  = *(p[k3]) && m3pre;
  st[10] = *(p[l3]) && m3pre;
  st[11] = *(p[m3]) && m3pre;
  st[12] = *(p[n3]) && m3pre;
  st[13] = *(p[o3]) && m3pre;
  st[14] = *(p[p3]) && n3pre;
  st[15] = *(p[q3]) && o3pre;
  st[16] = *(p[r3]) && p3pre;
  st[17] = *(p[s3]) && q3pre;
  st[18] = *(p[t3]) && r3pre;
  st[19] = *(p[u3]) && s3pre;
  st[20] = *(p[v3]) && t3pre;
  st[21] = *(p[w3]) && u3pre;
  st[22] = *(p[x3]) && v3pre;
  
  *showermax = 0;
  for (short i = 0; i < nShowerSums; i++) {
    short shower = (short) bsum(&st[i], nInShower);
    if (shower > *showermax)
      *showermax = shower;
  }

  // We need to find out if this gap was skipped.
  // The first step is to define variables which indicate
  // whether the set tubes on this gap which are projected 
  // to by each tube on the previous gap are MT.
  // =====================================================
  bool g2sl = *(p[f3]) || *(p[g3]) || *(p[h3]) ||
    (swi == 2 && (*(p[e3]) || *(p[i3]))) ||
    (swi == 3 && (*(p[d3]) || *(p[j3])));
         
  bool h2sl = *(p[g3]) || *(p[h3]) || *(p[i3]) ||
    (swi == 2 && (*(p[f3]) || *(p[j3]))) ||
    (swi == 3 && (*(p[e3]) || *(p[k3])));
         
  bool i2sl = *(p[h3]) || *(p[i3]) || *(p[j3]) ||
    (swi == 2 && (*(p[g3]) || *(p[k3]))) ||
    (swi == 3 && (*(p[f3]) || *(p[l3])));
         
  bool j2sl = *(p[i3]) || *(p[j3]) || *(p[k3]) ||
    (swi == 2 && (*(p[h3]) || *(p[l3]))) ||
    (swi == 3 && (*(p[g3]) || *(p[m3])));
         
  bool k2sl = *(p[j3]) || *(p[k3]) || *(p[l3]) ||
    (swi == 2 && (*(p[i3]) || *(p[m3]))) ||
    (swi == 3 && (*(p[h3]) || *(p[n3])));
         
  bool l2sl = *(p[k3]) || *(p[l3]) || *(p[m3]) ||
    (swi == 2 && (*(p[j3]) || *(p[n3]))) ||
    (swi == 3 && (*(p[i3]) || *(p[o3])));
         
  bool m2sl = *(p[l3]) || *(p[m3]) || *(p[n3]) ||
    (swi == 2 && (*(p[k3]) || *(p[o3]))) ||
    (swi == 3 && (*(p[j3]) || *(p[p3])));
         
  bool n2sl = *(p[m3]) || *(p[n3]) || *(p[o3]) ||
    (swi == 2 && (*(p[l3]) || *(p[p3]))) ||
    (swi == 3 && (*(p[k3]) || *(p[q3])));
         
  bool o2sl = *(p[n3]) || *(p[o3]) || *(p[p3]) ||
    (swi == 2 && (*(p[m3]) || *(p[q3]))) ||
    (swi == 3 && (*(p[l3]) || *(p[r3])));
         
  bool p2sl = *(p[o3]) || *(p[p3]) || *(p[q3]) ||
    (swi == 2 && (*(p[n3]) || *(p[r3]))) ||
    (swi == 3 && (*(p[m3]) || *(p[s3])));
         
  bool q2sl = *(p[p3]) || *(p[q3]) || *(p[r3]) ||
    (swi == 2 && (*(p[o3]) || *(p[s3]))) ||
    (swi == 3 && (*(p[n3]) || *(p[t3])));
         
  bool r2sl = *(p[q3]) || *(p[r3]) || *(p[s3]) ||
    (swi == 2 && (*(p[p3]) || *(p[t3]))) ||
    (swi == 3 && (*(p[o3]) || *(p[u3])));
         
  bool s2sl = *(p[r3]) || *(p[s3]) || *(p[t3]) ||
    (swi == 2 && (*(p[q3]) || *(p[u3]))) ||
    (swi == 3 && (*(p[p3]) || *(p[v3])));

  // If a tube on the current gap is pointed to by a tube 
  // the previous gap, then the current-gap tube's signal projects forward to 
  // the next gap iff all tubes pointed to by the previous-gap tube are MT.
  // ========================================================================
  bool d3g2skip3 = d3g2c && !g2sl;
  bool e3g2skip3 = e3g2c && !g2sl;
  bool f3g2skip3 = f3g2c && !g2sl;
  bool g3g2skip3 = g3g2c && !g2sl;
  bool h3g2skip3 = h3g2c && !g2sl;
  bool i3g2skip3 = i3g2c && !g2sl;
  bool j3g2skip3 = j3g2c && !g2sl;
                             
  bool e3h2skip3 = e3h2c && !h2sl;
  bool f3h2skip3 = f3h2c && !h2sl;
  bool g3h2skip3 = g3h2c && !h2sl;
  bool h3h2skip3 = h3h2c && !h2sl;
  bool i3h2skip3 = i3h2c && !h2sl;
  bool j3h2skip3 = j3h2c && !h2sl;
  bool k3h2skip3 = k3h2c && !h2sl;
                             
  bool f3i2skip3 = f3i2c && !i2sl;
  bool g3i2skip3 = g3i2c && !i2sl;
  bool h3i2skip3 = h3i2c && !i2sl;
  bool i3i2skip3 = i3i2c && !i2sl;
  bool j3i2skip3 = j3i2c && !i2sl;
  bool k3i2skip3 = k3i2c && !i2sl;
  bool l3i2skip3 = l3i2c && !i2sl;
                             
  bool g3j2skip3 = g3j2c && !j2sl;
  bool h3j2skip3 = h3j2c && !j2sl;
  bool i3j2skip3 = i3j2c && !j2sl;
  bool j3j2skip3 = j3j2c && !j2sl;
  bool k3j2skip3 = k3j2c && !j2sl;
  bool l3j2skip3 = l3j2c && !j2sl;
  bool m3j2skip3 = m3j2c && !j2sl;
                             
  bool h3k2skip3 = h3k2c && !k2sl;
  bool i3k2skip3 = i3k2c && !k2sl;
  bool j3k2skip3 = j3k2c && !k2sl;
  bool k3k2skip3 = k3k2c && !k2sl;
  bool l3k2skip3 = l3k2c && !k2sl;
  bool m3k2skip3 = m3k2c && !k2sl;
  bool n3k2skip3 = n3k2c && !k2sl;
                             
  bool i3l2skip3 = i3l2c && !l2sl;
  bool j3l2skip3 = j3l2c && !l2sl;
  bool k3l2skip3 = k3l2c && !l2sl;
  bool l3l2skip3 = l3l2c && !l2sl;
  bool m3l2skip3 = m3l2c && !l2sl;
  bool n3l2skip3 = n3l2c && !l2sl;
  bool o3l2skip3 = o3l2c && !l2sl;
                             
  bool j3m2skip3 = j3m2c && !m2sl;
  bool k3m2skip3 = k3m2c && !m2sl;
  bool l3m2skip3 = l3m2c && !m2sl;
  bool m3m2skip3 = m3m2c && !m2sl;
  bool n3m2skip3 = n3m2c && !m2sl;
  bool o3m2skip3 = o3m2c && !m2sl;
  bool p3m2skip3 = p3m2c && !m2sl;
                             
  bool k3n2skip3 = k3n2c && !n2sl;
  bool l3n2skip3 = l3n2c && !n2sl;
  bool m3n2skip3 = m3n2c && !n2sl;
  bool n3n2skip3 = n3n2c && !n2sl;
  bool o3n2skip3 = o3n2c && !n2sl;
  bool p3n2skip3 = p3n2c && !n2sl;
  bool q3n2skip3 = q3n2c && !n2sl;
                             
  bool l3o2skip3 = l3o2c && !o2sl;
  bool m3o2skip3 = m3o2c && !o2sl;
  bool n3o2skip3 = n3o2c && !o2sl;
  bool o3o2skip3 = o3o2c && !o2sl;
  bool p3o2skip3 = p3o2c && !o2sl;
  bool q3o2skip3 = q3o2c && !o2sl;
  bool r3o2skip3 = r3o2c && !o2sl;
                             
  bool m3p2skip3 = m3p2c && !p2sl;
  bool n3p2skip3 = n3p2c && !p2sl;
  bool o3p2skip3 = o3p2c && !p2sl;
  bool p3p2skip3 = p3p2c && !p2sl;
  bool q3p2skip3 = q3p2c && !p2sl;
  bool r3p2skip3 = r3p2c && !p2sl;
  bool s3p2skip3 = s3p2c && !p2sl;
                              
  bool n3q2skip3 = n3q2c && !q2sl;
  bool o3q2skip3 = o3q2c && !q2sl;
  bool p3q2skip3 = p3q2c && !q2sl;
  bool q3q2skip3 = q3q2c && !q2sl;
  bool r3q2skip3 = r3q2c && !q2sl;
  bool s3q2skip3 = s3q2c && !q2sl;
  bool t3q2skip3 = t3q2c && !q2sl;
                             
  bool o3r2skip3 = o3r2c && !r2sl;
  bool p3r2skip3 = p3r2c && !r2sl;
  bool q3r2skip3 = q3r2c && !r2sl;
  bool r3r2skip3 = r3r2c && !r2sl;
  bool s3r2skip3 = s3r2c && !r2sl;
  bool t3r2skip3 = t3r2c && !r2sl;
  bool u3r2skip3 = u3r2c && !r2sl;
                             
  bool p3s2skip3 = p3s2c && !s2sl;
  bool q3s2skip3 = q3s2c && !s2sl;
  bool r3s2skip3 = r3s2c && !s2sl;
  bool s3s2skip3 = s3s2c && !s2sl;
  bool t3s2skip3 = t3s2c && !s2sl;
  bool u3s2skip3 = u3s2c && !s2sl;
  bool v3s2skip3 = v3s2c && !s2sl;

  bool d3skip3 = d3g2skip3; 
  bool e3skip3 = e3g2skip3 || e3h2skip3;
  bool f3skip3 = f3g2skip3 || f3h2skip3 || f3i2skip3;
  bool g3skip3 = g3g2skip3 || g3h2skip3 || g3i2skip3 || g3j2skip3;
  bool h3skip3 = h3g2skip3 || h3h2skip3 || h3i2skip3 || h3j2skip3 || h3k2skip3;
  bool i3skip3 = i3g2skip3 || i3h2skip3 || i3i2skip3 || i3j2skip3 || i3k2skip3 || i3l2skip3;
  bool j3skip3 = j3g2skip3 || j3h2skip3 || j3i2skip3 || j3j2skip3 || j3k2skip3 || j3l2skip3 || j3m2skip3;
  bool k3skip3 = k3h2skip3 || k3i2skip3 || k3j2skip3 || k3k2skip3 || k3l2skip3 || k3m2skip3 || k3n2skip3;
  bool l3skip3 = l3i2skip3 || l3j2skip3 || l3k2skip3 || l3l2skip3 || l3m2skip3 || l3n2skip3 || l3o2skip3;
  bool m3skip3 = m3j2skip3 || m3k2skip3 || m3l2skip3 || m3m2skip3 || m3n2skip3 || m3o2skip3 || m3p2skip3;
  bool n3skip3 = n3k2skip3 || n3l2skip3 || n3m2skip3 || n3n2skip3 || n3o2skip3 || n3p2skip3 || n3q2skip3;
  bool o3skip3 = o3l2skip3 || o3m2skip3 || o3n2skip3 || o3o2skip3 || o3p2skip3 || o3q2skip3 || o3r2skip3;
  bool p3skip3 = p3m2skip3 || p3n2skip3 || p3o2skip3 || p3p2skip3 || p3q2skip3 || p3r2skip3 || p3s2skip3;
  bool q3skip3 = q3n2skip3 || q3o2skip3 || q3p2skip3 || q3q2skip3 || q3r2skip3 || q3s2skip3;
  bool r3skip3 = r3o2skip3 || r3p2skip3 || r3q2skip3 || r3r2skip3 || r3s2skip3;       
  bool s3skip3 = s3p2skip3 || s3q2skip3 || s3r2skip3 || s3s2skip3;
  bool t3skip3 = t3q2skip3 || t3r2skip3 || t3s2skip3;
  bool u3skip3 = u3r2skip3 || u3s2skip3;
  bool v3skip3 = v3s2skip3;
  
  // Projection to the next gap.
  // ===========================
  d_3 = d3c || d3skip3;
  e_3 = e3c || e3skip3;
  f_3 = f3c || f3skip3;
  g_3 = g3c || g3skip3;
  h_3 = h3c || h3skip3;
  i_3 = i3c || i3skip3;
  j_3 = j3c || j3skip3;
  k_3 = k3c || k3skip3;
  l_3 = l3c || l3skip3;
  m_3 = m3c || m3skip3;
  n_3 = n3c || n3skip3;
  o_3 = o3c || o3skip3;
  p_3 = p3c || p3skip3;
  q_3 = q3c || q3skip3;
  r_3 = r3c || r3skip3;
  s_3 = s3c || s3skip3;
  t_3 = t3c || t3skip3;
  u_3 = u3c || u3skip3;
  v_3 = v3c || v3skip3;
  
  return;
}



void 
mMuiFastRoadFinder::gap4(bool **p, 
                         short swi, short *showermax, bool *roadfound)
{

  // These variables indicate whether a tube in this gap is
  // pointed to by a tube in the previous gap. 
  // ======================================================
  bool a4d3 = d_3 && swi==3;
  bool b4d3 = d_3 && swi>=2;
  bool c4d3 = d_3;             
  bool d4d3 = d_3;             
  bool e4d3 = d_3;             
  bool f4d3 = d_3 && swi>=2;
  bool g4d3 = d_3 && swi==3;
                   
  bool b4e3 = e_3 && swi==3;
  bool c4e3 = e_3 && swi>=2;
  bool d4e3 = e_3;             
  bool e4e3 = e_3;             
  bool f4e3 = e_3;             
  bool g4e3 = e_3 && swi>=2;
  bool h4e3 = e_3 && swi==3;
                   
  bool c4f3 = f_3 && swi==3;
  bool d4f3 = f_3 && swi>=2;
  bool e4f3 = f_3;             
  bool f4f3 = f_3;             
  bool g4f3 = f_3;             
  bool h4f3 = f_3 && swi>=2;
  bool i4f3 = f_3 && swi==3;
                   
  bool d4g3 = g_3 && swi==3;
  bool e4g3 = g_3 && swi>=2;
  bool f4g3 = g_3;             
  bool g4g3 = g_3;             
  bool h4g3 = g_3;             
  bool i4g3 = g_3 && swi>=2;
  bool j4g3 = g_3 && swi==3;
                   
  bool e4h3 = h_3 && swi==3;
  bool f4h3 = h_3 && swi>=2;
  bool g4h3 = h_3;             
  bool h4h3 = h_3;             
  bool i4h3 = h_3;             
  bool j4h3 = h_3 && swi>=2;
  bool k4h3 = h_3 && swi==3;
                   
  bool f4i3 = i_3 && swi==3;
  bool g4i3 = i_3 && swi>=2;
  bool h4i3 = i_3;             
  bool i4i3 = i_3;             
  bool j4i3 = i_3;             
  bool k4i3 = i_3 && swi>=2;
  bool l4i3 = i_3 && swi==3;
                   
  bool g4j3 = j_3 && swi==3;
  bool h4j3 = j_3 && swi>=2;
  bool i4j3 = j_3;             
  bool j4j3 = j_3;             
  bool k4j3 = j_3;             
  bool l4j3 = j_3 && swi>=2;
  bool m4j3 = j_3 && swi==3;
                   
  bool h4k3 = k_3 && swi==3;
  bool i4k3 = k_3 && swi>=2;
  bool j4k3 = k_3;             
  bool k4k3 = k_3;             
  bool l4k3 = k_3;             
  bool m4k3 = k_3 && swi>=2;
  bool n4k3 = k_3 && swi==3;
                   
  bool i4l3 = l_3 && swi==3;
  bool j4l3 = l_3 && swi>=2;
  bool k4l3 = l_3;             
  bool l4l3 = l_3;             
  bool m4l3 = l_3;             
  bool n4l3 = l_3 && swi>=2;
  bool o4l3 = l_3 && swi==3;
                   
  bool j4m3 = m_3 && swi==3;
  bool k4m3 = m_3 && swi>=2;
  bool l4m3 = m_3;             
  bool m4m3 = m_3;             
  bool n4m3 = m_3;             
  bool o4m3 = m_3 && swi>=2;
  bool p4m3 = m_3 && swi==3;
                   
  bool k4n3 = n_3 && swi==3;
  bool l4n3 = n_3 && swi>=2;
  bool m4n3 = n_3;             
  bool n4n3 = n_3;             
  bool o4n3 = n_3;             
  bool p4n3 = n_3 && swi>=2;
  bool q4n3 = n_3 && swi==3;
                   
  bool l4o3 = o_3 && swi==3;
  bool m4o3 = o_3 && swi>=2;
  bool n4o3 = o_3;             
  bool o4o3 = o_3;             
  bool p4o3 = o_3;             
  bool q4o3 = o_3 && swi>=2;
  bool r4o3 = o_3 && swi==3;
                   
  bool m4p3 = p_3 && swi==3;
  bool n4p3 = p_3 && swi>=2;
  bool o4p3 = p_3;             
  bool p4p3 = p_3;             
  bool q4p3 = p_3;             
  bool r4p3 = p_3 && swi>=2;
  bool s4p3 = p_3 && swi==3;

  bool n4q3 = q_3 && swi==3;
  bool o4q3 = q_3 && swi>=2;
  bool p4q3 = q_3;             
  bool q4q3 = q_3;             
  bool r4q3 = q_3;             
  bool s4q3 = q_3 && swi>=2;
  bool t4q3 = q_3 && swi==3;
                   
  bool o4r3 = r_3 && swi==3;
  bool p4r3 = r_3 && swi>=2;
  bool q4r3 = r_3;             
  bool r4r3 = r_3;             
  bool s4r3 = r_3;             
  bool t4r3 = r_3 && swi>=2;
  bool u4r3 = r_3 && swi==3;
                   
  bool p4s3 = s_3 && swi==3;
  bool q4s3 = s_3 && swi>=2;
  bool r4s3 = s_3;             
  bool s4s3 = s_3;             
  bool t4s3 = s_3;             
  bool u4s3 = s_3 && swi>=2;
  bool v4s3 = s_3 && swi==3;

  bool q4t3 = t_3 && swi==3;
  bool r4t3 = t_3 && swi>=2;
  bool s4t3 = t_3;             
  bool t4t3 = t_3;             
  bool u4t3 = t_3;             
  bool v4t3 = t_3 && swi>=2;
  bool w4t3 = t_3 && swi==3;
                   
  bool r4u3 = u_3 && swi==3;
  bool s4u3 = u_3 && swi>=2;
  bool t4u3 = u_3;             
  bool u4u3 = u_3;             
  bool v4u3 = u_3;             
  bool w4u3 = u_3 && swi>=2;
  bool x4u3 = u_3 && swi==3;
                   
  bool s4v3 = v_3 && swi==3;
  bool t4v3 = v_3 && swi>=2;
  bool u4v3 = v_3;             
  bool v4v3 = v_3;             
  bool w4v3 = v_3;             
  bool x4v3 = v_3 && swi>=2;
  bool y4v3 = v_3 && swi==3;
                   
  bool a4pre = a4d3; 
  bool b4pre = b4d3 || b4e3;
  bool c4pre = c4d3 || c4e3 || c4f3;
  bool d4pre = d4d3 || d4e3 || d4f3 || d4g3;
  bool e4pre = e4d3 || e4e3 || e4f3 || e4g3 || e4h3;
  bool f4pre = f4d3 || f4e3 || f4f3 || f4g3 || f4h3 || f4i3;
  bool g4pre = g4d3 || g4e3 || g4f3 || g4g3 || g4h3 || g4i3 || g4j3;
  bool h4pre = h4e3 || h4f3 || h4g3 || h4h3 || h4i3 || h4j3 || h4k3;
  bool i4pre = i4f3 || i4g3 || i4h3 || i4i3 || i4j3 || i4k3 || i4l3;
  bool j4pre = j4g3 || j4h3 || j4i3 || j4j3 || j4k3 || j4l3 || j4m3;
  bool k4pre = k4h3 || k4i3 || k4j3 || k4k3 || k4l3 || k4m3 || k4n3;
  bool l4pre = l4i3 || l4j3 || l4k3 || l4l3 || l4m3 || l4n3 || l4o3;
  bool m4pre = m4j3 || m4k3 || m4l3 || m4m3 || m4n3 || m4o3 || m4p3;
  bool n4pre = n4k3 || n4l3 || n4m3 || n4n3 || n4o3 || n4p3 || n4q3;
  bool o4pre = o4l3 || o4m3 || o4n3 || o4o3 || o4p3 || o4q3 || o4r3;
  bool p4pre = p4m3 || p4n3 || p4o3 || p4p3 || p4q3 || p4r3 || p4s3;
  bool q4pre = q4n3 || q4o3 || q4p3 || q4q3 || q4r3 || q4s3 || q4t3;
  bool r4pre = r4o3 || r4p3 || r4q3 || r4r3 || r4s3 || r4t3 || r4u3;
  bool s4pre = s4p3 || s4q3 || s4r3 || s4s3 || s4t3 || s4u3 || s4v3;
  bool t4pre = t4q3 || t4r3 || t4s3 || t4t3 || t4u3 || t4v3;
  bool u4pre = u4r3 || u4s3 || u4t3 || u4u3 || u4v3;       
  bool v4pre = v4s3 || v4t3 || v4u3 || v4v3;
  bool w4pre = w4t3 || w4u3 || w4v3;
  bool x4pre = x4u3 || x4v3;
  bool y4pre = y4v3;
  
  // Search for road in current gap.
  // ===============================
  bool a4c = *(p[a4]) && a4pre;
  bool b4c = *(p[b4]) && b4pre;
  bool c4c = *(p[c4]) && c4pre;
  bool d4c = *(p[d4]) && d4pre;
  bool e4c = *(p[e4]) && e4pre;
  bool f4c = *(p[f4]) && f4pre;
  bool g4c = *(p[g4]) && g4pre;
  bool h4c = *(p[h4]) && h4pre;
  bool i4c = *(p[i4]) && i4pre;
  bool j4c = *(p[j4]) && j4pre;
  bool k4c = *(p[k4]) && k4pre;
  bool l4c = *(p[l4]) && l4pre;
  bool m4c = *(p[m4]) && m4pre;
  bool n4c = *(p[n4]) && n4pre;
  bool o4c = *(p[o4]) && o4pre;
  bool p4c = *(p[p4]) && p4pre;
  bool q4c = *(p[q4]) && q4pre;
  bool r4c = *(p[r4]) && r4pre;
  bool s4c = *(p[s4]) && s4pre;
  bool t4c = *(p[t4]) && t4pre;
  bool u4c = *(p[u4]) && u4pre;
  bool v4c = *(p[v4]) && v4pre;
  bool w4c = *(p[w4]) && w4pre;
  bool x4c = *(p[x4]) && x4pre;
  bool y4c = *(p[y4]) && y4pre;

  *roadfound = a4c || b4c || c4c || d4c || e4c || f4c || g4c || h4c || i4c || j4c || k4c || l4c || m4c 
                   || n4c || o4c || p4c || q4c || r4c || s4c || t4c || u4c || v4c || w4c || x4c || y4c;
  
  // Find shower size.
  // =================
  static const short nInShower = 5;
  static const short nShowerSums = 25;
  bool st[nShowerSums+nInShower-1];
  st[0]  = *(p[aaa4]) && a4pre;
  st[1]  = *(p[aa4])  && b4pre;
  st[2]  = *(p[a4])   && c4pre;
  st[3]  = *(p[b4])   && d4pre;
  st[4]  = *(p[c4])   && e4pre;
  st[5]  = *(p[d4])   && f4pre;
  st[6]  = *(p[e4])   && g4pre;
  st[7]  = *(p[f4])   && h4pre;
  st[8]  = *(p[g4])   && i4pre;
  st[9]  = *(p[h4])   && j4pre;
  st[10] = *(p[i4])   && k4pre;
  st[11] = *(p[j4])   && l4pre;
  st[12] = *(p[k4])   && m4pre;
  st[13] = *(p[l4])   && m4pre;
  st[14] = *(p[m4])   && m4pre;
  st[15] = *(p[n4])   && m4pre;
  st[16] = *(p[o4])   && m4pre;
  st[17] = *(p[p4])   && n4pre;
  st[18] = *(p[q4])   && o4pre;
  st[19] = *(p[r4])   && p4pre;
  st[20] = *(p[s4])   && q4pre;
  st[21] = *(p[t4])   && r4pre;
  st[22] = *(p[u4])   && s4pre;
  st[23] = *(p[v4])   && t4pre;
  st[24] = *(p[w4])   && u4pre;
  st[25] = *(p[x4])   && v4pre;
  st[26] = *(p[y4])   && w4pre;
  st[27] = *(p[z4])   && x4pre;
  st[28] = *(p[zz4])  && y4pre;
  
  *showermax = 0;
  for (short i = 0; i < nShowerSums; i++) {
    short shower = (short) bsum(&st[i], nInShower);
    if (shower > *showermax)
      *showermax = shower;
  }

  return;
}


short 
mMuiFastRoadFinder::bsum (bool *st, long n) 
{
  short sum = 0;
  for (short i = 0; i < n; i++) {
    if (st[i])
      sum++;
  }
  return sum;
}

---