AlgTrackCamList.cxx

Go to the documentation of this file.

// Package: CandTrackCam
//
// AlgTrackCamList.cxx
//
// marshall@hep.phy.cam.ac.uk
//
// chapman@hep.phy.cam.ac.uk
// blake@hep.phy.cam.ac.uk

#include <cassert>
#include <cmath>

#include "Conventions/Munits.h"
#include "MessageService/MsgService.h"
#include "MinosObjectMap/MomNavigator.h"
#include "Navigation/NavKey.h"
#include "Navigation/NavSet.h"
#include "Validity/VldContext.h"

#include "Algorithm/AlgFactory.h"
#include "Algorithm/AlgHandle.h"
#include "Algorithm/AlgConfig.h"
#include "Candidate/CandContext.h"

#include "RecoBase/CandTrackHandle.h"
#include "RecoBase/CandTrackListHandle.h"
#include "RecoBase/CandSliceHandle.h"
#include "RecoBase/CandSliceListHandle.h"
#include "RecoBase/CandStripHandle.h"
#include "RecoBase/CandStripListHandle.h"
#include "CandDigit/CandDeMuxDigitHandle.h"
#include "CandDigit/CandDeMuxDigitListHandle.h"

#include "CandTrackCam/AlgTrackCamList.h"
#include "CandTrackCam/CandTrackCamListHandle.h"
#include "CandTrackCam/CandTrackCamHandle.h"

#include "HitCam.h"
#include "ClusterCam.h"
#include "TrackSegmentCam.h"
#include "TrackCam.h"

#include "TObjArray.h"

ClassImp(AlgTrackCamList)

  CVSID("$Id: AlgTrackCamList.cxx,v 1.14 2007/11/11 08:48:25 rhatcher Exp $");

AlgTrackCamList::AlgTrackCamList() :
  StripWidth(4.108e-2)
{
}


AlgTrackCamList::~AlgTrackCamList()
{
}


void AlgTrackCamList::RunAlg(AlgConfig &ac,CandHandle &ch, CandContext &cx)
{
  assert(cx.GetDataIn());

  // Check for CandSliceListHandle input
  if (cx.GetDataIn()->InheritsFrom("CandSliceListHandle")) {
    const CandSliceListHandle *cslh = dynamic_cast<const CandSliceListHandle*>(cx.GetDataIn());
    const MomNavigator *mom = cx.GetMom();
    
    CandRecord* candrec = (CandRecord*)(cx.GetCandRecord());
    assert(candrec);
    
    vldc = (VldContext*)(candrec->GetVldContext());
    
    
    // Create the new tracklist
    CandTrackCamListHandle& tracklist = dynamic_cast<CandTrackCamListHandle&>(ch);
    
    if( !cslh || cslh->GetNDaughters()<1 ) {
      // Require number of CandSlices to be non-zero to do anything more
      MSG("AlgTrackCamList", Msg::kWarning) << " !cslh || cslh->GetNDaughters()<1 " << endl;
      return;
    }
 
    // Configure for correct detector
    switch(vldc->GetDetector()){
    case Detector::kFar:
      ModuleType=1; NumModules=2; PlanesInModule=248; break;
    case Detector::kNear:
      ModuleType=2; NumModules=1; PlanesInModule=120; break;
    case Detector::kCalDet:
      ModuleType=3; NumModules=1; PlanesInModule=60; break;
    default:
      ModuleType=-1; NumModules=0; PlanesInModule=0; break;
    }
     
 
    // Unpack AlgConfig
    CambridgeAnalysis=false;
    if( ac.KeyExists("CambridgeAnalysis") ) {CambridgeAnalysis = ac.GetInt("CambridgeAnalysis");}
  
 
    // Iterate over the slices
    TIter sliceItr(cslh->GetDaughterIterator());
     
    while (CandSliceHandle *slice = dynamic_cast<CandSliceHandle*>(sliceItr())) {

      // Find tracks in slice
      if(ModuleType==1 || ModuleType==3) {PECut=1.8; PECut2=1.8;}
      else if(ModuleType==2) {PECut=2.; PECut2=0.5;}
      RunTheFinder(slice);

        
      // Setup AlgTrackCam
      AlgFactory &af = AlgFactory::GetInstance();
      TString algtrackconfig = "default";
      if( ac.KeyExists("AlgTrackCamConfig") ) {
        const char* tmps;
        if( ac.Get("AlgTrackCamConfig",tmps)) { algtrackconfig = tmps; }
      }
      AlgHandle ah_trk = af.GetAlgHandle("AlgTrackCam", algtrackconfig.Data());

        
      // Create the CandContext
      CandContext cx0(this, mom);
      cx0.SetCandRecord(candrec);
        
      unsigned int NTracksFound = FinalTrackBank[0].size();

      // Run the finder again with different PECut settings
      if(CambridgeAnalysis==false && NTracksFound==0) {
        ClearUp(); PECut=0.5; PECut2=0.5;
        RunTheFinder(slice); 
        NTracksFound = FinalTrackBank[0].size();
      }
  
      // Sanity Check
      if(NTracksFound!= FinalTrackBank[1].size()) {
        MSG("AlgTrackCamList", Msg::kError) <<"Different Number of U and V Tracks Found!"<<endl;
        if( NTracksFound>=FinalTrackBank[1].size()) {NTracksFound = FinalTrackBank[1].size();}
      }

      //Loop over the U and V track found
      for(unsigned int i=0; i<NTracksFound; ++i) {
        TrackCam* tracku = FinalTrackBank[0][i];
        TrackCam* trackv = FinalTrackBank[1][i];
        
        //Fill the CandContext for this track
        TObjArray trackin;
        trackin.Add(tracku);
        trackin.Add(trackv);
        trackin.Add(slice);
        cx0.SetDataIn(&trackin);
          
        //Create the CandTrack
        CandTrackCamHandle cth = CandTrackCam::MakeCandidate(ah_trk, cx0);
        cth.SetName(TString("CandTrackCamHandle"));
        cth.SetTitle(TString("Created by AlgTrackCamList"));
          
        // Add CandTrack to CandTrackList
        if(cth.GetNDaughters()>0) {tracklist.AddDaughterLink(cth);}
      }
      
      ClearUp();

    } // End loop over slices
      
  } // End if inherits from CandSliceListHandle
  
}




void AlgTrackCamList::RunTheFinder(CandSliceHandle* slice)
{
  // Configure algorithm for the relevant detector
  // For ND, initial algorithm is applied only to calorimeter,
  // spectrometer is treated separately
  MSG("AlgTrackCamList", Msg::kVerbose) << "FINDING TRACKS IN SLICE" << endl;

  // Run the methods
  FormTheHits(slice);
  FormTheClusters(); 
  IDTrkAndShwClusters();
  FormTriplets();
  FindAllAssociations();
  FindPreferredJoins();
  FindMatchedJoins();
  FirstUVComparison();
  if(ModuleType==2) {NearDetectorTriplets();}
  Form2DTracks();
  Join2DTracks();
  SecondUVComparison();
  Form3DTracks(slice);

  return;
}




void AlgTrackCamList::FormTheHits(CandSliceHandle* slice)
{
  // Create and (store pointers to) HitCam objects for all strips in slice
  // meeting PH and XTalk requirements. HitCam class is essentially a wrapper
  // for CandStripHandles

  bool IsXTalk[2];
  int Counter; double StripCharge;
  
  // Loop over all strips in the slice
  TIter SliceItr(slice->GetDaughterIterator());
  while(CandStripHandle* SlcStrip = dynamic_cast<CandStripHandle*>(SliceItr())) {
    
    // Default Xtalk flag to false, as cannot identify ND XTalk
    IsXTalk[0]=false; IsXTalk[1]=false;
  
    // Identify XTalk strips for FD
    if(ModuleType==1) {
      IsXTalk[0]=true; IsXTalk[1]=true; Counter=0;
      
      TIter DigItr(SlcStrip->GetDaughterIterator());
      while(CandDeMuxDigitHandle* Digit = (CandDeMuxDigitHandle*)(DigItr())) {
        if( !((Digit->GetDeMuxDigitFlagWord()<8)
              && (Digit->GetDeMuxDigitFlagWord() & CandDeMuxDigit::kXTalk)==(CandDeMuxDigit::kXTalk)) ) {
          IsXTalk[Counter]=false;
        }
        if(Counter==0) {Counter=1;}
      }
    }
    
    // Add the hits to the bank, where they are stored in plane order
    if(IsXTalk[0]==false || IsXTalk[1]==false) {
      StripCharge=SlcStrip->GetCharge();

      if(StripCharge>PECut2 || StripCharge>PECut) {
        HitCam* hit = new HitCam(SlcStrip);

        if(StripCharge>PECut) {HitBank[SlcStrip->GetPlane()].push_back(hit);}
        if(StripCharge>PECut2) {AllHitBank[SlcStrip->GetPlane()].push_back(hit);}
      }
    }
  }
  

  // Print out the list of hits
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF HITS *** " << endl;
  for(int i=0; i<500; ++i){
    for(unsigned int j=0; j<HitBank[i].size(); ++j) {
      MSG("AlgTrackCamList", Msg::kVerbose)
        << " pln="  << HitBank[i][j]->GetPlane() 
        << " strp=" << HitBank[i][j]->GetStrip() 
        << " chg="  << HitBank[i][j]->GetCharge()
        << " time=" << HitBank[i][j]->GetTime() << endl;
    }
  }

  return;
}




void AlgTrackCamList::FormTheClusters()
{
  // For each plane where we stored hits, we look to form 1D clusters 
  // of these hits. This is controlled by the IsHitAssoc method in ClusterCam.
  // Pointers to the ClusterCam objects are stored in plane order in a 
  // clusterbank
  

  bool AddingHits;
  int i;

  for(int Module=0; Module<NumModules; ++Module) {
    for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
      
      // Loop over the hits on plane i
      for(unsigned int j=0; j<HitBank[i].size(); ++j) {
        
        // Make a cluster for the first hit on the plane
        if(HitBank[i][j]->GetUID()==0) {
          ClusterCam* Clust = new ClusterCam(HitBank[i][j]);
          ClusterBank[HitBank[i][j]->GetPlane()].push_back(Clust);

          // Change UID from 0 to 1 to show that hit has been added       
          HitBank[i][j]->SetUID(1); 
          AddingHits=true;
          
          // Loop over other hits on plane to form clusters
          while(AddingHits==true) {
            AddingHits=false;
            
            for(unsigned int k=0; k<HitBank[i].size(); ++k) {
              if(HitBank[i][k]->GetUID()==0 && Clust->IsHitAssoc(HitBank[i][k])) {
                Clust->AddHit(HitBank[i][k]); 
                HitBank[i][k]->SetUID(1);
                AddingHits=true;
              }
            }
          }

        }
      } // End loop over hits on plane
    } // End loop over planes in module
  } // End loop over modules

  return;
}




void AlgTrackCamList::IDTrkAndShwClusters()
{
  // Look at the 1D clusters formed and use simple techniques to get an idea
  // of how track-like or shower-like each cluster is.
  // Then look at all clusters on the plane and see how track-like or shower-like
  // the plane is.

  //  Detector::Detector_t Detector = vldc->GetDetector();

  int i, k0;
  int nclust0, nclust1, nhits0, nhits1, ShwAssocNum;
  vector<ClusterCam*> TempClust0;
  vector<ClusterCam*> TempClust1;


  // Simple checks to set initial track flags
  for(int Module=0; Module<NumModules; ++Module) {
    for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
      
      for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
        if(ClusterBank[i][j]->GetCharge()>PECut || ClusterBank[i][j]->GetDigits()>1 ) {
          ClusterBank[i][j]->SetTrkFlag(1);
        }
      }
    }
  }


  
  // Identify the shower-like clusters
  for(int Module=0; Module<NumModules; ++Module) {
    for(int Plane=1; Plane<PlanesInModule+1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
      
      for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
        ClusterCam* Clust0 = ClusterBank[i][j];

        nhits0=0; nhits1=0; nclust0=0; nclust1=0;
        TempClust0.clear(); TempClust1.clear(); 
        
        
        // Look for shower associations in nearby planes
        // Loop over nearby planes
        for(int k=-1; k<2; ++k) {
          k0=i+(2*k);
          
          // Check k0 is a valid nearby plane
          if(k0>Module*(PlanesInModule+1) && k0<(Module+1)*(PlanesInModule+1)) {

            // Loop over clusters on plane k0
            for(unsigned int k1=0; k1<ClusterBank[k0].size(); ++k1) {
              ClusterCam* Clust1 = ClusterBank[k0][k1];

              // Check shower-like associations
              if(Clust0->GetCharge()>1. && Clust1->GetCharge()>1.) {
                ShwAssocNum=Clust0->IsShwAssoc(Clust1);
                
                // Store the association data and the clusters
                if(ShwAssocNum>0) {
                  nhits0+=Clust1->GetEntries(); 
                  nclust0+=1;
                  TempClust0.push_back(Clust1);
                }
                
                if(ShwAssocNum>1) {
                  nhits1+=Clust1->GetEntries(); 
                  nclust1+=1;
                  TempClust1.push_back(Clust1);
                }
              }
            } // End loop over clusters on nearby plane
          }
        } // End loop over nearby planes


        // Use the above results to set the shower flags
        // Has important implications for finding tracks embedded in showers
        if(nclust1>1 && nhits1>4) {
          for(unsigned int k=0; k<TempClust1.size(); ++k) {
            TempClust1[k]->SetShwFlag(1);
          }
        }

        if(nclust0>4 && nhits0>5) {
          for(unsigned int k=0; k<TempClust0.size(); ++k) {
            TempClust0[k]->SetShwFlag(1);
          }
        }

      }
    }
  }



  // Identify track-like and shower-like planes
  double Charge;

  for(int Module=0; Module<NumModules; ++Module) {
    for(int Plane=0; Plane<PlanesInModule+1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
      Charge=0.;

      // Get total charge in clusters on plane
      for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
        Charge+=ClusterBank[i][j]->GetCharge();
      }

      if(Charge>0.) {
        for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
          ClusterCam* Clust = ClusterBank[i][j];

          // Set the plane flags for the clusters // should be <3
          if(Clust->GetEntries()<7 && (Clust->GetCharge()/Charge)>0.5) {Clust->SetTrkPlnFlag(1);}
          if(Clust->GetEntries()>1 && (Clust->GetCharge()/Charge)>0.5) {Clust->SetShwPlnFlag(1);}
        }
      }
    }
  }
  


  // Print out list of hits and 1D clusters
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF CLUSTERS *** " << endl;
  for(int i=0; i<500; ++i){
    for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
      MSG("AlgTrackCamList", Msg::kVerbose)
        << " plane="    << ClusterBank[i][j]->GetPlane() 
        << " begtpos="  << ClusterBank[i][j]->GetBegTPos() 
        << " endtpos="  << ClusterBank[i][j]->GetEndTPos() 
        << " trkflag="  << ClusterBank[i][j]->GetTrkFlag()
        << " shwpln="   << ClusterBank[i][j]->GetShwPlnFlag() 
        << " trkpln="   << ClusterBank[i][j]->GetTrkPlnFlag()
        << " shwflag="  << ClusterBank[i][j]->GetShwFlag() << endl;
    }
  }

  return;
}




void AlgTrackCamList::FormTriplets()
{
  // Treating U and V views separately, we consider associations between clusters
  // on nearby planes. We look for groups of three clusters, each on different 
  // planes, that can be joined together in a small track segment, or triplet. 

  // All possible triplets of the following form are created:
  
  // If we sit on plane 0, and there are no possible hits on a plane marked X, m 
  // means a lower plane number than our origin plane, and p means a higher plane 
  // number.
   
  // m 0 p         The simplest triplet
  
  // m X 0 p       One gap
  // m 0 X p
  
  // m X X 0 p     Two gaps
  // m X 0 X p
  // m 0 X X p


  int TripAssocNum;
  bool AlternateTriplets;
  bool JoinFlag;
  int i;
  
  // Begin loop to make the triplets
  for(int Module=0; Module<NumModules; ++Module) {

    for(int Plane=2; Plane<PlanesInModule-1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);

      // Loop over clusters on plane. This is our plane 0, the origin.
      for(unsigned int k0=0; k0<ClusterBank[i].size(); ++k0) {

        if(ClusterBank[i][k0]->GetTrkFlag()>0) {
          JoinFlag=false;


          // Look for triplets of form m 0 p
          if((i-2)>=0 && ClusterBank[i-2].size()>0 && (i+2)<500 && ClusterBank[i+2].size()>0) {
          
            // Look at the clusters on plane i-2...
            for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
            
              // ...and the clusters on plane i+2
              for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
              
                if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {

                  // Check the track-like association of the three clusters
                  TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][kp]);

                  // If the association is good, make the triplet
                  if(TripAssocNum>0) {
                    TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km], ClusterBank[i][k0], ClusterBank[i+2][kp]);

                    SegmentBank[i].push_back(seg0);
                    JoinFlag=true;

                    // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                    if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                  }
                }
              }
            }
          }
        

        
          // Look for triplets of form m <-> 0 p
          if((i+2)<500 && ClusterBank[i+2].size()>0) {

        
            // Look for triplets of form m X 0 p  
            if( Plane>3 && (i-4)>=0 && ClusterBank[i-4].size()>0) {

              // Look at the clusters on plane i+2...
              for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
              
                // ...and the clusters on plane i-4
                for(unsigned int km=0; km<ClusterBank[i-4].size(); ++km) {

                  if(ClusterBank[i-4][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {

                    // Check the track-like association of the three clusters
                    TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+2][kp]);
                  
                    // If the association is good, check whether we can make the triplet
                    if(TripAssocNum>0) {
                      AlternateTriplets=false;
                      
                      // Check to see if there is also an alternate triplet possibility with clusters on plane i-2
                      // i.e. Want to make sure that the X in "m X 0 p" is correct.

                      if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
                        // Look at the clusters on plane i-2
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][kp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      
                      // If everything is good, make the triplet
                      if(AlternateTriplets==false) {
                        TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-4][km],ClusterBank[i][k0],ClusterBank[i+2][kp]);
                        
                        SegmentBank[i].push_back(seg0);
                        JoinFlag=true;
        
                        // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                        if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                      }
                    }
                  }
                }
              }
            }


            // Look for triplets of form m X X 0 p  
            if( Plane>5 && (i-6)>=0 && ClusterBank[i-6].size()>0) {

              // Look at the clusters on plane i+2...
              for(unsigned int kp=0; kp<ClusterBank[i+2].size(); ++kp) {
              
                // ...and look at the clusters on plane i-6
                for(unsigned int km=0; km<ClusterBank[i-6].size(); ++km) {
                  
                  if(ClusterBank[i-6][km]->GetTrkFlag()>0 && ClusterBank[i+2][kp]->GetTrkFlag()>0) {

                    // Check the track-like association of the three clusters
                    TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i+2][kp]);
        
                    // If the association is good, check whether we can make the triplet
                    if(TripAssocNum>0) {
                      AlternateTriplets=false;
                      
                      // Check to see if there are alternate triplet possibilities with clusters on plane i-2 or i-4
                      // i.e. Want to make sure that the Xs in "m X X 0 p" are correct.

                      // Check first X
                      if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
                        // Look at any clusters on plane i-2
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i][k0])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][kp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      // Check second X, plane i-4
                      if(AlternateTriplets==false && ClusterBank[i-4].size()>0) {
                        // Look at any clusters on plane i-4
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i-4].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i-4][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i-4][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i][k0])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][ktmp],ClusterBank[i+2][kp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      // Check both Xs together, planes i-2 and i-4
                      if(AlternateTriplets==false && ClusterBank[i-2].size()>0 && ClusterBank[i-4].size()>0) {
                        // Look again at any clusters on plane i-4...
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i-4].size(); ++ktmp) {
                          
                          // ...and look again any clusters on plane i-2
                          for(unsigned int ktmp1=0; ktmp1<ClusterBank[i-2].size(); ++ktmp1) {
                            if(AlternateTriplets==false 
                               && ClusterBank[i-4][ktmp]->GetTrkFlag()>0 && ClusterBank[i-2][ktmp1]->GetTrkFlag()>0
                               && ClusterBank[i-4][ktmp]->IsTrkAssoc(ClusterBank[i-6][km],ClusterBank[i-2][ktmp1])>0
                               && ClusterBank[i-2][ktmp1]->IsTrkAssoc(ClusterBank[i-4][ktmp],ClusterBank[i][k0])>0
                               && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp1],ClusterBank[i+2][kp])>0)
                              {AlternateTriplets=true;}
                          }
                          
                        }
                      }
                      
                      // If everything is good, make the triplet
                      if(AlternateTriplets==false) {
                        TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-6][km],ClusterBank[i][k0],ClusterBank[i+2][kp]);
                        
                        SegmentBank[i].push_back(seg0);
                        JoinFlag=true;

                        // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                        if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                      }
                    }
                  }
                }
              }
            }

          }



          // Look for triplets of form m 0 <-> p
          if((i-2)>0 && ClusterBank[i-2].size()>0) {

        
            // Look for triplets of form m 0 X p  
            if( Plane<PlanesInModule-3 && (i+4)<500 && ClusterBank[i+4].size()>0) {

              // Look at the clusters on plane i-2...
              for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
              
                // ... and the clusters on plane i+4
                for(unsigned int kp=0; kp<ClusterBank[i+4].size(); ++kp) {

                  if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+4][kp]->GetTrkFlag()>0) {

                    // Check the track-like association of the three clusters
                    TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+4][kp]);
                   
                    // If the association is good, check whether we can make the triplet
                    if(TripAssocNum>0) {
                      AlternateTriplets=false;
                      
                      // Check to see if there is also an alternate triplet possibility with clusters on plane i+2
                      // i.e. Want to make sure that the X in "m 0 X p" is correct.

                      if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
                        // Look at the clusters on plane i+2
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      
                      // If everything is good, make the triplet
                      if(AlternateTriplets==false) {
                        // Store segment on empty plane too
                        for(int j=0; j<2; ++j) {
                          TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km],ClusterBank[i][k0],ClusterBank[i+4][kp]);
                          
                          SegmentBank[i+2*j].push_back(seg0);
                          JoinFlag=true;
                        }
                        
                        // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                        if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                      }
                    }
                  }
                }
              }
            }


            // Look for triplets of form m 0 X X p  
            if(Plane<PlanesInModule-5 && (i+6)<500 && ClusterBank[i+6].size()>0) {

              // Look at the clusters on plane i-2...
              for(unsigned int km=0; km<ClusterBank[i-2].size(); ++km) {
              
                // ...and look at the clusters on plane i+6
                for(unsigned int kp=0; kp<ClusterBank[i+6].size(); ++kp) {

                  if(ClusterBank[i-2][km]->GetTrkFlag()>0 && ClusterBank[i+6][kp]->GetTrkFlag()>0) {
                    // Check the track-like association of the three clusters
                    TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+6][kp]);
                    
                    // If the association is good, check whether we can make the triplet
                    if(TripAssocNum>0) {
                      AlternateTriplets=false;
        
                      // Check to see if there are alternate triplet possibilities with clusters on plane i+2 or i+4
                      // i.e. Want to make sure that the Xs in "m 0 X X p" are correct.
              
                      // Check first X, plane i+2
                      if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
                        // Look at any clusters on plane i+2
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+6][kp])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      // Check second X, plane i+4
                      if(AlternateTriplets==false && ClusterBank[i+4].size()>0) {
                        // Look at any clusters on plane i+4
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i+4].size(); ++ktmp) {
                          if(AlternateTriplets==false && ClusterBank[i+4][ktmp]->GetTrkFlag()>0
                             && ClusterBank[i+4][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+6][kp])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+4][ktmp])>0) 
                            {AlternateTriplets=true;}
                        }
                      }
                      
                      // Check both Xs together, planes i+2 and i+4
                      if(AlternateTriplets==false && ClusterBank[i+2].size()>0 && ClusterBank[i+4].size()>0) {
                        // Look again at any clusters on plane i+2...
                        for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
                          
                          // ...and look again at any clusters on plane i+4
                          for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+4].size(); ++ktmp1) {
                            if(AlternateTriplets==false 
                               && ClusterBank[i+2][ktmp]->GetTrkFlag()>0 &&ClusterBank[i+4][ktmp1]->GetTrkFlag()>0
                               && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][km],ClusterBank[i+2][ktmp])>0
                               && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][ktmp1])>0
                               && ClusterBank[i+4][ktmp1]->IsTrkAssoc(ClusterBank[i+2][ktmp],ClusterBank[i+6][kp])>0)
                              {AlternateTriplets=true;}
                          }
                          
                        }
                      }
                      
                      // If everything is good, make the triplet
                      if(AlternateTriplets==false) {
                        // Store segment on empty planes too
                        for(int j=0; j<3; ++j) {
                          TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-2][km],ClusterBank[i][k0],ClusterBank[i+6][kp]);
                          
                          SegmentBank[i+2*j].push_back(seg0);
                          JoinFlag=true;
                        }
                        // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                        if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                      }
                    }
                  }
                }
              }
            }

          }



          // Look for triplets of form m X 0 X p
          if(Plane>3 && Plane<PlanesInModule-3 && (i-4)>=0 && ClusterBank[i-4].size()>0 && (i+4)<500 && ClusterBank[i+4].size()>0) {

            // Look at clusters on planes i-4 and i+4
            for(unsigned int km=0; km<ClusterBank[i-4].size(); ++km) {
              for(unsigned int kp=0; kp<ClusterBank[i+4].size(); ++kp) {

                if(ClusterBank[i-4][km]->GetTrkFlag()>0 && ClusterBank[i+4][kp]->GetTrkFlag()>0) {
                  // Check the track-like association of the three clusters
                  TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+4][kp]);
                  
                  // If the association is good, check whether we can make the triplet
                  if(TripAssocNum>0) {
                    AlternateTriplets=false;
                    
                    // Check to see if there are alternate triplet possibilities with clusters on plane i-2 or i+2
                    // i.e. Want to make sure that the Xs in "m X 0 X p" are correct.
                    
                    // Check first X, plane i-2
                    if(AlternateTriplets==false && ClusterBank[i-2].size()>0) {
                      // Look at any clusters on plane i-2
                      for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
                        if(AlternateTriplets==false && ClusterBank[i-2][ktmp]->GetTrkFlag()>0
                           && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
                           && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+4][kp])>0) 
                          {AlternateTriplets=true;}
                      }
                    }
                    
                    // Check second X, plane i+2
                    if(AlternateTriplets==false && ClusterBank[i+2].size()>0) {
                      // Look at any clusters on plane i+2
                      for(unsigned int ktmp=0; ktmp<ClusterBank[i+2].size(); ++ktmp) {
                        if(AlternateTriplets==false  && ClusterBank[i+2][ktmp]->GetTrkFlag()>0
                           && ClusterBank[i+2][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0 
                           && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i+2][ktmp])>0) 
                          {AlternateTriplets=true;}
                      }
                    }
                    
                    // Check both Xs together, planes i-2 and i+2
                    if(AlternateTriplets==false && ClusterBank[i-2].size()>0 && ClusterBank[i+2].size()>0) {
                      
                      // Loop again at any clusters on plane i-2
                      for(unsigned int ktmp=0; ktmp<ClusterBank[i-2].size(); ++ktmp) {
                        
                        // Look again at any clusters on plane i+2
                        for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+2].size(); ++ktmp1) {
                          if(AlternateTriplets==false 
                             && ClusterBank[i-2][ktmp]->GetTrkFlag()>0 && ClusterBank[i+2][ktmp1]->GetTrkFlag()>0
                             && ClusterBank[i-2][ktmp]->IsTrkAssoc(ClusterBank[i-4][km],ClusterBank[i][k0])>0 
                             && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-2][ktmp],ClusterBank[i+2][ktmp1])>0 
                             && ClusterBank[i+2][ktmp1]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+4][kp])>0)
                            {AlternateTriplets=true;}
                        }
                      }
                    }

                    // If everything is good, make the triplet
                    if(AlternateTriplets==false) {
                      // Store segment on empty planes too
                      for(int j=0; j<2; ++j) {
                        TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-4][km],ClusterBank[i][k0],ClusterBank[i+4][kp]);
                        
                        SegmentBank[i+2*j].push_back(seg0);
                        JoinFlag=true;
                      }
                      // Indicate that triplet is exceptionally track-like, by setting flag to 2.
                      if(TripAssocNum>1) {ClusterBank[i][k0]->SetTrkFlag(2);}
                    }               
                  }
                }
              }
            }
          }


          // Set the exceptionally track-like flag for lone hits that form part of a triplet
          if(JoinFlag==true && (ClusterBank[i][k0]->GetEndTPos()==ClusterBank[i][k0]->GetBegTPos())) {
            ClusterBank[i][k0]->SetTrkFlag(2);
          }

          
        }

      } // End loop over clusters on plane
    
    } // End loop over planes in module

  } // End loop over modules


  // Print out list of hits and 1D clusters
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** 2ND LIST OF CLUSTERS *** " << endl;
  for(int i=0; i<500; ++i){
    for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
      MSG("AlgTrackCamList", Msg::kVerbose)
        << " plane="    << ClusterBank[i][j]->GetPlane() 
        << " begtpos="  << ClusterBank[i][j]->GetBegTPos() 
        << " endtpos="  << ClusterBank[i][j]->GetEndTPos() 
        << " trkflag="  << ClusterBank[i][j]->GetTrkFlag()
        << " shwpln="   << ClusterBank[i][j]->GetShwPlnFlag() 
        << " trkpln="   << ClusterBank[i][j]->GetTrkPlnFlag()
        << " shwflag="  << ClusterBank[i][j]->GetShwFlag() << endl;
    }
  }


  // Print out list of triplets
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF TRIPLETS *** " << endl;
  for(int i=0; i<500; ++i) {
    for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
      MSG("AlgTrackCamList", Msg::kVerbose) << " Segment Num " << j << ", Plane " << i <<endl;
 
      ClusterCam* clr0 = SegmentBank[i][j]->GetCluster(0);
      ClusterCam* clr1 = SegmentBank[i][j]->GetCluster(1);
      ClusterCam* clr2 = SegmentBank[i][j]->GetCluster(2);

      MSG("AlgTrackCamList", Msg::kVerbose) 
        << " (" << clr0->GetPlane() << "," << clr0->GetBegTPos() << "->" << clr0->GetEndTPos() << ") "
        << " (" << clr1->GetPlane() << "," << clr1->GetBegTPos() << "->" << clr1->GetEndTPos() << ") "
        << " (" << clr2->GetPlane() << "," << clr2->GetBegTPos() << "->" << clr2->GetEndTPos() << ") " << endl;
    }
  }
  
  return;
}




void AlgTrackCamList::FindAllAssociations()
{
  // For each triplet formed, we consider nearby triplets and see whether
  // the triplets are associated. 

  // For each triplet, we simply find the nearby other triplets that have 
  // a compatible beginning/end position/direction.

  // Later we look more that one triplet away, to find strings of 
  // "preferred associations" that might indicate a track.


  int i=0;


  // Loop over the planes
  for(int Module=0; Module<NumModules; ++Module) {

    for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
        
      if((i+2)<500 && SegmentBank[i].size()>0 && SegmentBank[i+2].size()>0) {
        
        // Look at the segments on plane i
        for(unsigned int k0=0; k0<SegmentBank[i].size(); ++k0) {
          
          // For each of these, look at the segments on plane i+2
          for(unsigned int kp=0; kp<SegmentBank[i+2].size(); ++kp) {
            
            // Check associations. 
            if(SegmentBank[i][k0]->IsAssoc(SegmentBank[i+2][kp])) {
              
              // If segments are associated, we add to list of ALL possible joins
              SegmentBank[i][k0]->AddAssocSegToEnd(SegmentBank[i+2][kp]);
              SegmentBank[i+2][kp]->AddAssocSegToBeg(SegmentBank[i][k0]);
            }
          }
        }
      }
    } // End loop over planes in module
  } // End loop over modules




  // Print out list of all associations
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ALL SEG ASSOCIATIONS *** " << endl;

  // Print out list of triplets
  for(int i=0; i<500; ++i) {
    for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
      TrackSegmentCam* segment = SegmentBank[i][j];

      MSG("AlgTrackCamList", Msg::kVerbose) << "--- Plane " << i << ", Seg Number " << j 
                                            << ", BegPlane " << segment->GetBegPlane() 
                                            << ", EndPlane " << segment->GetEndPlane() 
                                            << ", BegTPos " << segment->GetBegTPos() 
                                            << ", EndTPos " << segment->GetEndTPos() 
                                            <<endl;
 
      MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " <<endl;
      for(unsigned int k=0; k<segment->GetNAssocSegBeg(); ++k) {
        TrackSegmentCam* segtmp = segment->GetAssocSegBeg(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;
      }
      
      MSG("AlgTrackCamList", Msg::kVerbose) << " End: " <<endl;
      for(unsigned int k=0; k<segment->GetNAssocSegEnd(); ++k) {
        TrackSegmentCam* segtmp = segment->GetAssocSegEnd(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;           
      }
    }
  }

  return;
}




void AlgTrackCamList::FindPreferredJoins()
{
  // Having made all possible associations between triplets above, we 
  // select those associations that are most track-like and so are preferred.
  
  // For a given triplet, we know which triplets are associated with the
  // beginning and which are associated at the end. If there are associations
  // between these beginning and end triplets too, then we are quite likely 
  // to be considering track-like segments.



  int i;

  vector<TrackSegmentCam*> mTempSeg;
  vector<TrackSegmentCam*> pTempSeg;
  
  bool JoinFlag; bool mJnFlag; bool pJnFlag; 

  double GradientTolerance=8.;
  if(ModuleType==2) {GradientTolerance=5.;}


  // Loop over the planes
  for(int Module=0; Module<NumModules; ++Module) {

    for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
      i=Plane + Module*(PlanesInModule+1);
      

      // SET THE TEMPORARY TRACK FLAGS
     
      // Consider previous adjacent triplets   ( i-2 <-> ) i <-> i+2
      //                                                  Seg
      // If there are associations at both ends of triplet, set triplet's 
      // temporary track flag to 1. 

      // If one of beginning associated triplets also starts at lower plane 
      // than current triplet, set the triplet's temporary track flag to 2
      mTempSeg.clear();

      // Loop over the triplets centered on plane i
      for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
        TrackSegmentCam* Seg = SegmentBank[i][j];

        // Check to see if there are associations at end
        if(Seg->GetNAssocSegEnd()>0) {
          mTempSeg.push_back(Seg);
          
          // Check to see if there are associations at beginning
          if(Seg->GetNAssocSegBeg()>0) {
            
            // We now know that triplet has possible associations at both ends
            Seg->SetTmpTrkFlag(1);

            // Loop over the associations at beginning and set temp track flags
            // for segments on plane i
            for(unsigned int k=0; k<Seg->GetNAssocSegBeg(); ++k) {
              if(Seg->GetTmpTrkFlag()<2
                 && (Seg->GetAssocSegBeg(k)->GetBegPlane() < Seg->GetBegPlane()) ) 
                {Seg->SetTmpTrkFlag(2); break;}
            }
          }
        }
      }



      // Subsequent adjacent triplets   i <-> i+2 ( <-> i+4 )
      //                                      Segp
      // If there are associations at both ends of triplet centered on 
      // plane i+2, set triplet's temporary track flag to 1.

      // If one of the end associated triplets also ends at higher plane
      // than i+2 triplet, set i+2 triplet's temporary track flag to 2
      pTempSeg.clear();

      // Loop over the triplets centered on plane i+2
      for(unsigned int j=0; j<SegmentBank[i+2].size(); ++j) {
        TrackSegmentCam* Segp = SegmentBank[i+2][j];

        // Check to see if there are associations at beginning
        if(Segp->GetNAssocSegBeg()>0) {
          pTempSeg.push_back(Segp);
          
          // Check to see if there are associations at end
          if(Segp->GetNAssocSegEnd()>0) {
            
            // We now know that triplet has possible associations at both ends
            Segp->SetTmpTrkFlag(1);
            
            // Loop over associations at end and set temp track flags
            // for segments on plane i+2
            for(unsigned int k=0; k<Segp->GetNAssocSegEnd(); ++k) {
              if(Segp->GetTmpTrkFlag()<2
                 && (Segp->GetAssocSegEnd(k)->GetEndPlane() > Segp->GetEndPlane()) ) 
                {Segp->SetTmpTrkFlag(2); break;}
            }
          }
        }
      }
      
  

      // Look for preferred joins   ( i-2 <-> ) i <-> i+2
      // Loop the over segments on plane i+2, which we know have possible beginning associations
      for(unsigned int j=0; j<pTempSeg.size(); ++j) {
        JoinFlag=false;

        // Loop over these beginning associations. If one beginning associated
        // triplet already has temp track flag 2, set JoinFlag to true.
        for(unsigned int k=0; k<pTempSeg[j]->GetNAssocSegBeg(); ++k) {
          if(pTempSeg[j]->GetAssocSegBeg(k)->GetTmpTrkFlag()>1) {JoinFlag=true; break;}   
        }

        // If JoinFlag is true, and a beginning assoc segment has temp track flag 1 (i.e. possible 
        // beginning and end associations), set the temp track flag to 2 for this beg assoc segment
        if(JoinFlag==true) {
          // Loop over beginning associations again and set temp track flags
          // for segments on plane i
          for(unsigned int k=0; k<pTempSeg[j]->GetNAssocSegBeg(); ++k) {
            if(pTempSeg[j]->GetAssocSegBeg(k)->GetTmpTrkFlag()==1) {
              pTempSeg[j]->GetAssocSegBeg(k)->SetTmpTrkFlag(2); 
            }
          }
        }
      }



      // Look for preferred joins    i <-> i+2 ( <-> i+4 )
      // Loop over segments on plane i, which we know have possible end associations
      for(unsigned int j=0; j<mTempSeg.size(); ++j) {
        JoinFlag=false;

        // Loop over these end associations. If one end associated triplet already
        // has temp track flag 2, set JoinFlag to true.
        for(unsigned int k=0; k<mTempSeg[j]->GetNAssocSegEnd(); ++k) {
          if(mTempSeg[j]->GetAssocSegEnd(k)->GetTmpTrkFlag()>1) {JoinFlag=true; break;}   
        }
        
        // If JoinFlag is true, and an end assoc segment has temp track flag 1 (i.e. possible 
        // beginning and end associations), set the temp track flag to 2 for this end assoc segment
        if(JoinFlag==true) {
          // Loop over end associations again and set temp track flags
          // for segments on plane i+2
          for(unsigned int k=0; k<mTempSeg[j]->GetNAssocSegEnd(); ++k) {
            if(mTempSeg[j]->GetAssocSegEnd(k)->GetTmpTrkFlag()==1) {
              mTempSeg[j]->GetAssocSegEnd(k)->SetTmpTrkFlag(2); 
            }
          }
        }
      }


      // Now have temporary track flags set to 2 for all segments on planes i and i+2
      // which look like very promising track segments      




      // MAKE THE PREFERRED JOINS BETWEEN TRIPLETS ON PLANE i AND PLANE i+2
      if(mTempSeg.size()>0 && pTempSeg.size()>0) {
      
        // Look for doubly preferred joins    ( i-2 <-> ) i <-> i+2 ( <-> i+4 )

        // Loop over segments on plane i, which we know have possible end associations
        for(unsigned int j=0; j<mTempSeg.size(); ++j) {
          TrackSegmentCam* Seg = mTempSeg[j];

          // For these, loop over possible end associations
          for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
            TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
 
            // If both segments given temp track flag 2 above check assocations 
            // between triplets on either side
            if(Seg->GetTmpTrkFlag()>1 && Segp->GetTmpTrkFlag()>1) {

              // For each Segp, see if there is a Segm association
              mJnFlag=false;
              for(unsigned int l=0; l<Seg->GetNAssocSegBeg(); ++l) {
                TrackSegmentCam* Segm = Seg->GetAssocSegBeg(l);
                if(Segm->IsAssoc(Segp) ) {mJnFlag=true; break;}
              }

              // For each end assoc of Segp, see if it is assoc with Seg
              pJnFlag=false;
              for(unsigned int l=0; l<Segp->GetNAssocSegEnd(); ++l) {
                TrackSegmentCam* Segp2 = Segp->GetAssocSegEnd(l);
                if(Seg->IsAssoc(Segp2) ) {pJnFlag=true; break;}
              }

              // Make the preferred assocations
              if(mJnFlag==true && pJnFlag==true){
                if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
                        (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
                   && Seg->GetBegPlane()<=Segp->GetBegPlane()
                   && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
                  Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}

                else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
                                                            << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
                                                            << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
                                                            << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
              }
            }
          }
        }
     


        // Look for singly preferred joins (1)    ( i-2 <-> ) i <-> i+2 

        // Loop over segments on plane i, which we know have possible end associations
        for(unsigned int j=0; j<mTempSeg.size(); ++j) {
          TrackSegmentCam* Seg = mTempSeg[j];
          
          // If segment given temp track flag 2 above 
          if(Seg->GetTmpTrkFlag()>1) {
            JoinFlag=false;

            // Check to see if we will already have made the preferred association above
            for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
              TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
              
              if(Segp->GetTmpTrkFlag()>1) {JoinFlag=true; break;}
            }

            // If have not made the association, can check assocations 
            // between triplets on either side
            if(JoinFlag==false) {
              for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
                TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);

                // For each Segp, see if there is a Segm associated
                mJnFlag=false;
                for(unsigned int l=0; l<Seg->GetNAssocSegBeg(); ++l) {
                  TrackSegmentCam* Segm = Seg->GetAssocSegBeg(l);
                  if(mJnFlag==false && Segm->IsAssoc(Segp) ) {mJnFlag=true; break;}
                }
               
                // Make the preferred assocations
                if(mJnFlag==true){
                  if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
                          (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
                     && Seg->GetBegPlane()<=Segp->GetBegPlane()
                     && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
                    Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}

                  else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
                                                              << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
                                                              << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
                                                              << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
                }
              }
            }
          }
        }



        // Look for singly preferred joins (2)    i <-> i+2 ( <-> i+4 )

        // Loop the over segments on plane i+2, which we know have possible beginning associations
        for(unsigned int j=0; j<pTempSeg.size(); ++j) {
          TrackSegmentCam* Segp = pTempSeg[j];

          // If segment given temp track flag 2 above
          if(Segp->GetTmpTrkFlag()>1) {
            JoinFlag=false;

            // Check to see if we will already have made the preferred association above
            for(unsigned int k=0; k<Segp->GetNAssocSegBeg(); ++k) {
              TrackSegmentCam* Seg = Segp->GetAssocSegBeg(k);
              
              if(Seg->GetTmpTrkFlag()>1) {JoinFlag=true; break;}
            }

            // If have not made the association, can check assocations 
            // between triplets on either side
            if(JoinFlag==false) {
              for(unsigned int k=0; k<Segp->GetNAssocSegBeg(); ++k) {
                TrackSegmentCam* Seg = Segp->GetAssocSegBeg(k);

                // For each Seg, see if there is a Segp2 associated
                pJnFlag=false;
                for(unsigned int l=0; l<Segp->GetNAssocSegEnd(); ++l) {
                  TrackSegmentCam* Segp2 = Segp->GetAssocSegEnd(l);
                  if(pJnFlag==false && Seg->IsAssoc(Segp2) ) {pJnFlag=true; break;}
                }
               
                // Make the preferred assocations
                if(pJnFlag==true){
                  if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
                          (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
                     && Seg->GetBegPlane()<=Segp->GetBegPlane()
                     && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
                  Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}

                  else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
                                                              << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
                                                              << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
                                                              << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
                }
              }
            }
          }
        }



        // Look for other joins we have missed    i <-> i+2 

        // Loop over segments on plane i, which we know have possible end associations
        for(unsigned int j=0; j<mTempSeg.size(); ++j) {
          TrackSegmentCam* Seg = mTempSeg[j];
          
          // For these, loop over the end assocations
          for(unsigned int k=0; k<Seg->GetNAssocSegEnd(); ++k) {
            TrackSegmentCam* Segp = Seg->GetAssocSegEnd(k);
            
            // If segments and its end assoc segment were given temp track flag 2 above
            if(Seg->GetTmpTrkFlag()<2 && Segp->GetTmpTrkFlag()<2) {
              
              // Look at each possible end assoc for Seg to see if we've already
              // made the assocation
              mJnFlag=false;
              for(unsigned int l=0; l<Seg->GetNAssocSegEnd(); ++l) {
                if(Seg->GetAssocSegEnd(l)->GetTmpTrkFlag()>1) {mJnFlag=true; break;}
              }

              // Look at each possible beginning assoc for Segp to see if we've already
              // made the assocation
              pJnFlag=false;
              for(unsigned int l=0; l<Segp->GetNAssocSegBeg(); ++l) {
                if(Segp->GetAssocSegBeg(l)->GetTmpTrkFlag()>1) {pJnFlag=true; break;}
              }
                      
              // Make the preferred associations if it hasn't been done by one of the loops above
              if(mJnFlag==false && pJnFlag==false) {
                if(fabs((Seg->GetBegTPos()-Seg->GetEndTPos())-
                        (Segp->GetBegTPos()-Segp->GetEndTPos()))<GradientTolerance*StripWidth
                   && Seg->GetBegPlane()<=Segp->GetBegPlane()
                   && Seg->GetEndPlane()<=Segp->GetEndPlane() ) {
                Seg->AddPrefSegToEnd(Segp); Segp->AddPrefSegToBeg(Seg);}

                else {MSG("AlgTrackCamList", Msg::kVerbose) << "Seg " << Seg->GetBegPlane() << " " << Seg->GetEndPlane() << " "
                                                            << Seg->GetBegTPos() << " " << Seg->GetEndTPos() << " Segp " 
                                                            << Segp->GetBegPlane() << " " << Segp->GetEndPlane() << " "
                                                            << Segp->GetBegTPos() << " " << Segp->GetEndTPos() << endl;}
              }
            }       
          }
        }
        


      }



      // Clear flags ready to consider next planes
      for(unsigned int j=0; j<mTempSeg.size(); ++j) {
        mTempSeg[j]->SetTmpTrkFlag(0);
      }
      
      for(unsigned int j=0; j<pTempSeg.size(); ++j) {
        pTempSeg[j]->SetTmpTrkFlag(0);
      }


    } // End loop over planes in module

  } // End loop over modules

  
  
  
  // Look for any other preferred associations we may have missed.
  // Find segments with no preferred end association and look a bit
  // further to find preferred associations.

  // Not for ND
  if(ModuleType==1 || ModuleType==3) {

    int NewPlane; int Increment;
    double SegTPos, NearbySegTPos; 
    bool AssocsStopHere;
    bool ConsiderSegment;

    for(int Module=0; Module<NumModules; ++Module) {
      for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
        i=Plane + Module*(PlanesInModule+1);
          
        // Loop over all segments
        for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
          TrackSegmentCam* Seg0 = SegmentBank[i][j];
        
          // Loop over the segments with no preferred end association
          if(Seg0->GetNPrefSegEnd()==0 && Seg0->GetNPrefSegBeg()>0) {
          
            // Look at nearby segments to see if these already continue the associations.
            AssocsStopHere=true;
            SegTPos=0.5*(Seg0->GetBegTPos()+Seg0->GetEndTPos());

            for(unsigned int k=0; k<SegmentBank[i].size(); ++k) {
              TrackSegmentCam* NearbySeg = SegmentBank[i][k];
            
              if(NearbySeg==Seg0) {continue;}
            
              NearbySegTPos=0.5*(NearbySeg->GetBegTPos()+NearbySeg->GetEndTPos());

              if(fabs(NearbySegTPos-SegTPos)<1 && 
                 NearbySeg->GetNPrefSegEnd()>0 && NearbySeg->GetNPrefSegBeg()>0) {AssocsStopHere=false; break;}
            }

            if(AssocsStopHere==false) {break;}


            JoinFlag=false;

            // Consider the segment, plus all the segments, SegBeg, with which 
            // it has a preferred beginning association.
            for(int k=-1; k<int(Seg0->GetNPrefSegBeg()); ++k) {
              // First time, consider Seg0, rather than one of its preferred beginning associations.
              TrackSegmentCam* SegBeg = 0;
              if(k<0) {SegBeg=Seg0;}
              else{SegBeg = Seg0->GetPrefSegBeg(k);}
              Increment=4;

              // Look up to 12 planes ahead within the same module.
              while(Increment<=14) {
                NewPlane=SegBeg->GetEndPlane()+Increment;
              
                if(NewPlane>=(Module+1)*(PlanesInModule+1)) {break;}

                // See if SegBeg is associated with a segment on this plane, NextSeg.
                // If so, also see if SegBeg is associated with one of NextSeg's 
                // preferred end associated segments.

                // If so, make the preferred associations between SegBeg and NextSeg.
                for(unsigned int l=0; l<SegmentBank[NewPlane].size(); ++l) {
                  TrackSegmentCam* NextSeg = SegmentBank[NewPlane][l];


                  // Check segments are not now linked in a chain of preferred associations.
                  ConsiderSegment=true;
                  for(unsigned int m=0; m<NextSeg->GetNPrefSegBeg(); ++m) {
                    TrackSegmentCam* PrefSeg = NextSeg->GetPrefSegBeg(m);

                    if(PrefSeg->GetTmpTrkFlag()==1) {NextSeg->SetTmpTrkFlag(1); ConsiderSegment=false; break;}
                  }
                  if(ConsiderSegment==false) {continue;}


                  if( SegBeg->IsAssoc(NextSeg) ) {
                  
                    for(unsigned int m=0; m<NextSeg->GetNPrefSegEnd(); ++m) {
                      TrackSegmentCam* SegEnd = NextSeg->GetPrefSegEnd(m);
                    
                      if( SegBeg->IsAssoc(SegEnd) ) {
                        JoinFlag=true;
                      
                        SegBeg->AddPrefSegToEnd(NextSeg);
                        NextSeg->AddPrefSegToBeg(SegBeg);

                        SegBeg->SetTmpTrkFlag(1); NextSeg->SetTmpTrkFlag(1);

                        MSG("AlgTrackCamList", Msg::kVerbose) << "Made missing Pref assoc, SegBeg "
                                                              << SegBeg->GetBegPlane() << ", " << SegBeg->GetEndPlane() << ", "
                                                              << SegBeg->GetBegTPos() << ", " << SegBeg->GetEndTPos()
                                                              << " NextSeg " << NextSeg->GetBegPlane() << ", " 
                                                              << NextSeg->GetEndPlane() << ", " << NextSeg->GetBegTPos() 
                                                              << ", " << NextSeg->GetEndTPos() << endl;
                        break;
                      }
                    }
                  }
                }
                Increment+=2;
              }

              // Reset TmpTrkFlags
              for(int p=i; p<(1+Module)*PlanesInModule; ++p) {
                for(unsigned int q=0; q<SegmentBank[p].size(); ++q) {SegmentBank[p][q]->SetTmpTrkFlag(0);}
              }

              if(JoinFlag==true) {break;}
            }
          }
        

        
          // Now loop over the segments with no preferred beginning association
          if(Seg0->GetNPrefSegBeg()==0 && Seg0->GetNPrefSegEnd()>0) {
          
            // Look at nearby segments to see if these already continue the associations.
            AssocsStopHere=true;
            SegTPos=0.5*(Seg0->GetBegTPos()+Seg0->GetEndTPos());

            for(unsigned int k=0; k<SegmentBank[i].size(); ++k) {
              TrackSegmentCam* NearbySeg = SegmentBank[i][k];
            
              if(NearbySeg==Seg0) {continue;}
            
              NearbySegTPos=0.5*(NearbySeg->GetBegTPos()+NearbySeg->GetEndTPos());

              if(fabs(NearbySegTPos-SegTPos)<1 && 
                 NearbySeg->GetNPrefSegEnd()>0 && NearbySeg->GetNPrefSegBeg()>0) {AssocsStopHere=false; break;}
            }

            if(AssocsStopHere==false) {break;}


            JoinFlag=false;

            // Consider the segment, plus all the segments, SegEnd, with which 
            // it has a preferred end association.
            for(int k=-1; k<int(Seg0->GetNPrefSegEnd()); ++k) {
              // First time, consider Seg0, rather than one of its preferred end associations.
              TrackSegmentCam* SegEnd = 0;
              if(k<0) {SegEnd=Seg0;}
              else{SegEnd = Seg0->GetPrefSegEnd(k);}
              Increment=4;

              // Look up to 12 planes behind within the same module.
              while(Increment<=14) {
                NewPlane=SegEnd->GetBegPlane()-Increment;
              
                if(NewPlane<=(Module)*(PlanesInModule+1)) {break;}

                // See if SegEnd is associated with a segment on this plane, PrevSeg.
                // If so, also see if SegEnd is associated with one of PrevSeg's 
                // preferred beginning associated segments.

                // If so, make the preferred associations between SegEnd and PrevSeg.
                for(unsigned int l=0; l<SegmentBank[NewPlane].size(); ++l) {
                  TrackSegmentCam* PrevSeg = SegmentBank[NewPlane][l];


                  // Check segments are not now linked in a chain of preferred associations.
                  ConsiderSegment=true;
                  for(unsigned int m=0; m<PrevSeg->GetNPrefSegEnd(); ++m) {
                    TrackSegmentCam* PrefSeg = PrevSeg->GetPrefSegEnd(m);

                    if(PrefSeg->GetTmpTrkFlag()==1) {PrevSeg->SetTmpTrkFlag(1); ConsiderSegment=false; break;}
                  }
                  if(ConsiderSegment==false) {continue;}
        

                  if( PrevSeg->IsAssoc(SegEnd) ) {
                  
                    for(unsigned int m=0; m<PrevSeg->GetNPrefSegBeg(); ++m) {
                      TrackSegmentCam* SegBeg = PrevSeg->GetPrefSegBeg(m);
                            
                      if( SegBeg->IsAssoc(SegEnd) ) {
                        JoinFlag=true;
                      
                        SegEnd->AddPrefSegToBeg(PrevSeg);
                        PrevSeg->AddPrefSegToEnd(SegEnd);

                        SegEnd->SetTmpTrkFlag(1); PrevSeg->SetTmpTrkFlag(1);

                        MSG("AlgTrackCamList", Msg::kVerbose) << "Made missing Pref assoc, SegEnd "
                                                              << SegEnd->GetBegPlane() << ", " << SegEnd->GetEndPlane() << ", "
                                                              << SegEnd->GetBegTPos() << ", " << SegEnd->GetEndTPos()
                                                              << " PrevSeg " << PrevSeg->GetBegPlane() << ", " 
                                                              << PrevSeg->GetEndPlane() << ", " << PrevSeg->GetBegTPos() 
                                                              << ", " << PrevSeg->GetEndTPos() << endl;
                        break;
                      }
                    }
                  }
                }
                Increment+=2;
              }

              // Reset TmpTrkFlags
              for(int p=i; p>Module*(PlanesInModule); --p) {
                for(unsigned int q=0; q<SegmentBank[p].size(); ++q) {SegmentBank[p][q]->SetTmpTrkFlag(0);}
              }

              if(JoinFlag==true) {break;}
            }
          }

        }
      }
    }

  }
 
  


  // Print out list of preferred associations
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF PREF SEG ASSOCIATIONS *** " << endl;

  // Print out list of triplets
  for(int i=0; i<500; ++i) {
    for(unsigned int j=0;j<SegmentBank[i].size(); ++j) {
      TrackSegmentCam* segment = SegmentBank[i][j];

      MSG("AlgTrackCamList", Msg::kVerbose) << " ----Plane " << i << ", Seg Number " << j 
                                            << ", BegPlane " << segment->GetBegPlane() 
                                            << ", EndPlane " << segment->GetEndPlane() 
                                            << ", BegTPos " << segment->GetBegTPos() 
                                            << ", EndTPos " << segment->GetEndTPos() 
                                            <<endl;
 
      MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " << endl;
      for(unsigned int k=0; k<segment->GetNPrefSegBeg(); ++k) {
        TrackSegmentCam* segtmp = segment->GetPrefSegBeg(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Pref number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;
      }
      
      MSG("AlgTrackCamList", Msg::kVerbose) << " End: " << endl;
      for(unsigned int k=0; k<segment->GetNPrefSegEnd(); ++k) {
        TrackSegmentCam* segtmp = segment->GetPrefSegEnd(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Pref number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;           
      }
    }
  }
  
  return;
}




void AlgTrackCamList::FindMatchedJoins()
{
  // Having made all the preferred assocations, we actually match triplets 
  // together to form longer segments of the track.

  // We use the idea of seed segments to identify the last segment in a chain
  // of segments.



  int i;
  bool JoinFlag;

  // Loop over the planes
  for(int Module=0; Module<NumModules; ++Module) {
    for(int Plane=3; Plane<PlanesInModule-1; ++Plane) {
  
      i=Plane + Module*(PlanesInModule+1);
      
      // Loop over the triplets centered on plane i
      for(unsigned int j=0; j<SegmentBank[i].size(); ++j) {
        
        // Firstly, consider isolated segments. These must be their own seed segments.
        if(SegmentBank[i][j]->GetNAssocSegBeg()==0 && SegmentBank[i][j]->GetNAssocSegEnd()==0) {
          ViewSegBank[SegmentBank[i][j]->GetPlaneView()].push_back(SegmentBank[i][j]);
          SegmentBank[i][j]->SetSeedSegment(SegmentBank[i][j]);
        }


        // Then, consider those segments which have some preferred associations
        if(SegmentBank[i][j]->GetNPrefSegBeg()>0 || SegmentBank[i][j]->GetNPrefSegEnd()>0) {
          JoinFlag=false;

          // If there is one pref assoc at beginning (always false for first segment in module)
          if(SegmentBank[i][j]->GetNPrefSegBeg()==1) {
            TrackSegmentCam* Segm = SegmentBank[i][j]->GetPrefSegBeg(0);

            // Check that we are considering a simple chain of segments
            if(Segm->GetNPrefSegEnd()==1) {
              // Get the last segment in the chain of segments
              TrackSegmentCam* SeedSeg = Segm->GetSeedSegment();

              SeedSeg->AddSegment(SegmentBank[i][j]); 
              SegmentBank[i][j]->SetSeedSegment(SeedSeg);
              SegmentBank[i][j]->SetUID(-1);
              JoinFlag=true;
            }
          }

          // Now consider those triplets which aren't just one pref assoc at beg and one at end
          if(JoinFlag==false) {
            ViewSegBank[SegmentBank[i][j]->GetPlaneView()].push_back(SegmentBank[i][j]);
            SegmentBank[i][j]->SetSeedSegment(SegmentBank[i][j]);
            
            // Loop over the preferred associations at beginning
            for(unsigned int k=0; k<SegmentBank[i][j]->GetNPrefSegBeg(); ++k) {
              TrackSegmentCam* SeedSeg = SegmentBank[i][j]->GetPrefSegBeg(k)->GetSeedSegment();
              TrackSegmentCam* MatchSeg = SegmentBank[i][j];

              if(SeedSeg->GetBegPlane()<=MatchSeg->GetBegPlane() && SeedSeg->GetEndPlane()<=MatchSeg->GetEndPlane()) {

                if( (SeedSeg->GetEntries()>3 && MatchSeg->GetEntries()>3)
                    || 4.*fabs( (SeedSeg->GetBegTPos()-SeedSeg->GetEndTPos())/(SeedSeg->GetEndPlane()-SeedSeg->GetBegPlane()) -
                                (MatchSeg->GetBegTPos()-MatchSeg->GetEndTPos())/(MatchSeg->GetEndPlane()-MatchSeg->GetBegPlane()) ) < 10*StripWidth ) {
                  MatchSeg->AddMatchSegToBeg(SeedSeg);
                  SeedSeg->AddMatchSegToEnd(MatchSeg);
                }
              }

            }
          }
        }
      } // End loop over triplets centered on plane i
      
    } // End loop over planes in module

  } // End loop over modules



  // Form matched associations between non-adjacent segments
  vector<TrackSegmentCam*> TempSeg1;
  vector<TrackSegmentCam*> TempSeg2;

  const int npts0=4; int npts; double inv_npts; int Module; int Plane;
  bool CheckCoilHole;
  
  
  // Loop over the planeviews U and V
  for(int view=0; view<2; ++view) {
    TempSeg1.clear(); TempSeg2.clear();

 
    // Loop over entries stored in ViewSegBank, above
    for(unsigned int i=0; i<ViewSegBank[view].size(); ++i) {
      TrackSegmentCam* Seg1 = ViewSegBank[view][i];

      // If no matched entries at end of segment
      if(Seg1->GetNMatchSegEnd()==0) {
        JoinFlag=false;

        
        // For FD, attempt to match segments across the coil hole
        // For ND, need to be careful making associations, so we don't join two separate tracks together

        // Check which module we are in
        Module=int(Seg1->GetEndPlane()/(PlanesInModule+1));
        npts=npts0;

        // These are the tpos values required to locate the FD coil hole
        if( (ModuleType==1 && Seg1->GetEndTPos()>-0.62 && Seg1->GetEndTPos()<0.62) 
            || ( 1+(Seg1->GetEndPlane()-Seg1->GetBegPlane())/2>5) ) {
          
          // Calculate range of possible association and configure
          // - worst case scenario - track crosses entire width of coil hole, tcoil = ~0.6m
          // - it will take a distance in z, zcoil = tcoil * dz/dt to do this
          // - zcoil is equivalent to nplanes = zcoil/0.0594 = tcoil / ( 0.0594 * dt/dz )
          // - npts = nplanes/2.0 = 1.0 / ( 0.198 * dt/dz)
          // - where dt/dz = Seg1->GetEndDir();
          inv_npts=0.198*Seg1->GetEndDir(); //was inv_npts=0.5*Seg1->GetEndDir();
          
          if(inv_npts<0) {inv_npts=-inv_npts;}
          if(inv_npts<0.01) {inv_npts=0.01;}
          npts=int(1./inv_npts);

          if(npts<13) {npts=13;}    //replaced 8 with 13
          if(npts>100) {npts=100;}  //replaced 18 with 30 - possibly this should be a function of the gradient of the track?
        }

        
        // Loop over this range of possible association
        for(int j=1; j<npts; ++j) {
          Plane=Seg1->GetEndPlane()+(2*j);
          
          if( JoinFlag==false && Plane<((PlanesInModule+1)*(Module+1)) ) {
            
            // Loop over segments on each possible plane
            for(unsigned int k=0; k<SegmentBank[Plane].size(); ++k) {
              TrackSegmentCam* Seg2 = SegmentBank[Plane][k];

              // Check association, first across coil-hole 
              if(ModuleType==1 && (( Seg2->GetBegTPos()>-0.62 && Seg2->GetBegTPos()<0.21
                                     && Seg1->GetEndTPos()>-0.21 && Seg1->GetEndTPos()<0.62)
                                   || (Seg2->GetBegTPos()>-0.21 && Seg2->GetBegTPos()<0.62
                                       && Seg1->GetEndTPos()>-0.62 && Seg1->GetEndTPos()<0.21)) ) {CheckCoilHole=true;}
              else {CheckCoilHole=false;}
              

              // Then main association checks
              if( (Seg2->GetSeedSegment()==Seg2 && Seg2->GetNMatchSegBeg()==0) 
                  && (Seg2->GetBegPlane()-Seg1->GetEndPlane()>0 || (Seg1->GetEntries()>5 && Seg2->GetEntries()>5) )
                  && (j<npts0 || CheckCoilHole==true || (1+(Seg1->GetEndPlane()-Seg1->GetBegPlane())/2>5) ) ) {
                
                if( Seg1->IsAssoc(Seg2) ) {

                  // ND nearby 2D track protection
                  if( ModuleType!=2 || (fabs(Seg1->GetEndDir()-Seg2->GetBegDir())<0.25 && 
                                        ((Seg1->GetEntries()<11 && Seg2->GetEntries()<11 
                                          && (ModuleType!=2 || Seg2->GetBegPlane()-Seg1->GetEndPlane()<50)) 
                                         || Seg2->GetBegPlane()-Seg1->GetEndPlane()<11) ) ) {

                    TempSeg1.push_back(Seg1);
                    TempSeg2.push_back(Seg2);
                    
                    JoinFlag=true;
                    
                    MSG("AlgTrackCamList", Msg::kVerbose) << "   NON-ADJ ASSOC: " 
                                                          << Seg1->GetBegPlane() << " " << Seg1->GetEndPlane() << " (" 
                                                          << Seg1->GetBegTPos() << ", " << Seg1->GetEndTPos() << ") "
                                                          << " enddir " << Seg1->GetEndDir()
                                                          << " Entries " << Seg1->GetEntries() << " "
                                                          << Seg2->GetBegPlane() << " " << Seg2->GetEndPlane() << " (" 
                                                          << Seg2->GetBegTPos() << ", " << Seg2->GetEndTPos() << ") " 
                                                          << " begdir " << Seg2->GetBegDir() << " Entries " << Seg2->GetEntries() << endl; 
                  }
                }
              }
            }
          }
        } 
      }
    } // End loop over entries in ViewSegBank
    

    // Make the associations
    for(unsigned int i=0; i<TempSeg1.size(); ++i) {
      TrackSegmentCam* Seg1 = TempSeg1[i];
      TrackSegmentCam* Seg2 = TempSeg2[i];

      if(Seg1->GetBegPlane()<=Seg2->GetBegPlane() && Seg1->GetEndPlane()<=Seg2->GetEndPlane()) {
        
        if( (Seg1->GetEntries()>3 && Seg2->GetEntries()>3)
            || 4.*fabs( (Seg1->GetBegTPos()-Seg1->GetEndTPos())/(Seg1->GetBegPlane()-Seg1->GetEndPlane()) -
                        (Seg2->GetBegTPos()-Seg2->GetEndTPos())/(Seg2->GetBegPlane()-Seg2->GetEndPlane()) ) < 10*StripWidth ) {
          Seg2->AddMatchSegToBeg(Seg1);
          Seg1->AddMatchSegToEnd(Seg2);
        }
      }

    }

  } // End loop over planeviews



  // Print out list of segments
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF MATCHED SEG ASSOCIATIONS *** " << endl;

  for(int View=0; View<2; ++View) {

    MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;

    for(unsigned int i=0; i<ViewSegBank[View].size(); ++i) {
      TrackSegmentCam* segment = (TrackSegmentCam*)(ViewSegBank[View][i]);
      MSG("AlgTrackCamList", Msg::kVerbose) << "---  Seg0 number =" << i
                                            << "  begpln: "  << segment->GetBegPlane()
                                            << "  begtpos: " << segment->GetBegTPos()
                                            << "  endpln: "  << segment->GetEndPlane()
                                            << "  endtpos: " << segment->GetEndTPos() << endl;

      MSG("AlgTrackCamList", Msg::kVerbose) << "     Beg: " << endl;
      for(unsigned int j=0; j<segment->GetNMatchSegBeg(); ++j) {
        TrackSegmentCam* segtmp = segment->GetMatchSegBeg(j);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Match number=" << j
                                              << "  begpln: "  << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: "  << segtmp->GetEndPlane()
                                              << "  endtpos: "  << segtmp->GetEndTPos() << endl;
      }
      
      MSG("AlgTrackCamList", Msg::kVerbose) << "     End: " << endl;
      for(unsigned int j=0; j<segment->GetNMatchSegEnd(); ++j) {
        TrackSegmentCam* segtmp = segment->GetMatchSegEnd(j);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Match number=" << j
                                              << "  begpln: "  << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: "  << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;           
      }
    }
  }


}





void AlgTrackCamList::FirstUVComparison()
{
  // Having made sections of the track, from 'matched' segments, we compare
  // the sections in the U view with those in the V view, looking for compatibity.

  // First all clusters are investigated to mark the segments as track-like or
  // shower-like. Then track-like segments which have 'overlapping' partners in the 
  // opposite view are marked by setting their UID values to 2

  MSG("AlgTrackCamList", Msg::kVerbose) << " *** Comparing views (1) *** " << endl;


  int shwctr, plnctr;
  unsigned int nsegments=0;
  unsigned int nsegmentsV=0;
  vector<TrackSegmentCam*> mysegbnk[2];


  // Determine track-like segments
  for(int View=0; View<2; ++View) {

    // Loop over all segments in the current view
    nsegments = ViewSegBank[View].size();
    for(unsigned int i=0; i<nsegments; ++i) {

      TrackSegmentCam* Seg = ViewSegBank[View][i];
      shwctr=0; plnctr=0;
      
      // Loop over all the clusters in the current segment
      const unsigned int ncluster = Seg->GetEntries();
      for(unsigned int j=0;j<ncluster;++j) {
        ClusterCam* cluster = Seg->GetCluster(j);
        
        // Find how track-like or shower-like the cluster is
        if(cluster->GetTrkFlag()>1) { 
          if(cluster->GetShwFlag()<1) {shwctr++;}
        }
        if(cluster->GetTrkPlnFlag()>0) {plnctr++;}
      }
      
      // Set the UID accordingly. This depends greatly on the initial values 
      // of the shower flags used when identifying track and shower clusters.
      if(CambridgeAnalysis==false) {Seg->SetUID(1);}

      if( plnctr>1 && Seg->GetEntries()>4 ) { // Just added
        if(CambridgeAnalysis==true && shwctr>0) {Seg->SetUID(1);} 
        if(shwctr>2) {Seg->SetUID(2);} 
      }
    }
  }


  // Loop over segments in U view
  nsegments = ViewSegBank[0].size();
  for(unsigned int i=0; i<nsegments; ++i) {

    TrackSegmentCam* Segu = ViewSegBank[0][i];
    
    // For these, loop over segments in V view
    nsegmentsV = ViewSegBank[1].size();
    for(unsigned int j=0; j<nsegmentsV; ++j) {

      TrackSegmentCam* Segv = ViewSegBank[1][j];

      // Select pairs of segments, one from each view, which overlap by more than 1 plane
      // and which both have UID>0 (i.e. were marked as track-like above)
      if( Segu->GetBegPlane()+1<Segv->GetEndPlane() 
          && Segu->GetEndPlane()>Segv->GetBegPlane()+1 
          && Segu->GetUID()>0 && Segv->GetUID()>0 ) {
        if(Segv->GetUID()>1) {mysegbnk[0].push_back(Segu); }
        if(Segu->GetUID()>1) {mysegbnk[1].push_back(Segv); }
      }
    }
  }

  
  // Loop over the segments we just selected and set their UIDs to 2
  for(int View=0; View<2; ++View) {
    nsegments = mysegbnk[View].size();
    for(unsigned int i=0; i<nsegments; ++i) {
      mysegbnk[View][i]->SetUID(2);
    }
  }



  // Mark clusters that aren't going to be available for ND +/-10 plane track finding
  if(vldc->GetDetector()==Detector::kNear) {
    
    // First mark all those for which we will probably find a +/-2 plane track
    for(int View=0; View<2; ++View) {
      for(unsigned int i=0; i<mysegbnk[View].size(); ++i) {
        TrackSegmentCam* Seg = mysegbnk[View][i];

        for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
          Seg->GetCluster(j)->SetNDFlag(0);  
        }
      }
    }


    // Then look specifically for clusters that are in +/-2 plane sections of detector
    double BegTPos, EndTPos, Centre, HitTPos;
    double Tolerance=2.*StripWidth;

    // Loop over all clusters
    for(int i=0; i<PlanesInModule; ++i) {
      for(unsigned int j=0; j<ClusterBank[i].size(); ++j) {
        ClusterCam* Clust = ClusterBank[i][j];

        // Get cluster properties       
        BegTPos=Clust->GetBegTPos();
        EndTPos=Clust->GetEndTPos();
        Centre=0.5*(BegTPos+EndTPos);


        // Now compare with the hits +/-2 planes away
        if(Clust->GetNDFlag()>0 && i>1) {
          for(unsigned int k=0; k<HitBank[i-2].size(); ++k) {
            HitTPos=HitBank[i-2][k]->GetTPos();

            if(fabs(HitTPos-Centre)<Tolerance || fabs(HitTPos-BegTPos)<Tolerance 
               || fabs(HitTPos-EndTPos)<Tolerance) {Clust->SetNDFlag(0); break;}
          }
        }

        if(Clust->GetNDFlag()>0) {
          for(unsigned int k=0; k<HitBank[i+2].size(); ++k) {
            HitTPos=HitBank[i+2][k]->GetTPos();
            
            if(fabs(HitTPos-Centre)<Tolerance || fabs(HitTPos-BegTPos)<Tolerance 
               || fabs(HitTPos-EndTPos)<Tolerance) {Clust->SetNDFlag(0); break;}
          }
        }
      }
    }

  } // End ND only section


  // Print out the results
  for(int View=0; View<2; ++View) {
    MSG("AlgTrackCamList", Msg::kVerbose) << "  View =" << View << endl;

    nsegments = ViewSegBank[View].size();
    for(unsigned int i=0; i<nsegments; ++i) {

      TrackSegmentCam* Seg = ViewSegBank[View][i];

      MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg number=" << i 
                                            << "  begpln: "  << Seg->GetBegPlane()
                                            << "  begtpos: " << Seg->GetBegTPos()
                                            << "  endpln: "  << Seg->GetEndPlane()
                                            << "  endtpos: " << Seg->GetEndTPos() 
                                            << "  UID: "     << Seg->GetUID() 
                                            << "  clusters: "<< Seg->GetEntries()
                                            << endl;
    }
  }
  
  return;
}




void AlgTrackCamList::NearDetectorTriplets()
{
  MSG("AlgTrackCamList", Msg::kVerbose) << " MAKE ND TRIPLETS " << endl;

  int TripAssocNum; bool AlternateTriplets;
  vector<TrackSegmentCam*> NDViewSegBank[2];
  vector<TrackSegmentCam*> TempNDViewSegBank[2];


  // MAKE THE TRIPLETS
  for(int i=10; i<PlanesInModule-9; ++i) {

    // Only consider Fully Instrumented planes
    if((i-1)%10!=0 && (i-6)%10!=0) {continue;}

    // Loop over clusters on plane. This is our plane 0, the origin.
    for(unsigned int k0=0; k0<ClusterBank[i].size(); ++k0) {
    
      if(ClusterBank[i][k0]->GetTrkFlag()>0 && ClusterBank[i][k0]->GetNDFlag()>0) {
        
        // Look for triplets of form m 0 p
        if( (i-10)>=0 && ClusterBank[i-10].size()>0 && (i+10)<=PlanesInModule && ClusterBank[i+10].size()>0 ) {
          
          // Look at the clusters available on plane i-10...
          for(unsigned int km=0; km<ClusterBank[i-10].size(); ++km) {
            
            // ...and the clusters available on plane i+10
            for(unsigned int kp=0; kp<ClusterBank[i+10].size(); ++kp) {
              
              if( ClusterBank[i-10][km]->GetTrkFlag()>0 && ClusterBank[i-10][km]->GetNDFlag()>0
                  && ClusterBank[i+10][kp]->GetTrkFlag()>0 && ClusterBank[i+10][kp]->GetNDFlag()>0 ) {
                
                // Check the track-like association of the three clusters
                TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+10][kp],10);
        
                // If the association is good, make the triplet
                if(TripAssocNum>0) {
                  TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-10][km], ClusterBank[i][k0], ClusterBank[i+10][kp]);
                  for(unsigned int j=0; j<seg0->GetEntries(); ++j) {seg0->GetCluster(j)->SetNDFlag(2);}

                  NDSegmentBank[i].push_back(seg0);
                }
              }
            }
          }
        }



        
        // Look for triplets of form m X 0 p  
        if( (i-20)>=0 && ClusterBank[i-20].size()>0 && (i+10)<=PlanesInModule && ClusterBank[i+10].size()>0 ) {

          // Look at the clusters on plane i+10...
          for(unsigned int kp=0; kp<ClusterBank[i+10].size(); ++kp) {
              
            // ...and the clusters on plane i-20
            for(unsigned int km=0; km<ClusterBank[i-20].size(); ++km) {

              if( ClusterBank[i-20][km]->GetTrkFlag()>0 && ClusterBank[i+10][kp]->GetTrkFlag()>0
                  && ClusterBank[i-20][km]->GetNDFlag()>0 && ClusterBank[i+10][kp]->GetNDFlag()>0 ) {

                // Check the track-like association of the three clusters
                TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+10][kp],10);
                  
                // If the association is good, check whether we can make the triplet
                if(TripAssocNum>0) {
                  AlternateTriplets=false;
                      
                  // Check to see if there is also an alternate triplet possibility with clusters on plane i-2
                  // i.e. Want to make sure that the X in "m X 0 p" is correct.

                  if(AlternateTriplets==false && ClusterBank[i-10].size()>0) {
                    // Look at the clusters on plane i-10
                    for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
                      if(AlternateTriplets==false && ClusterBank[i-10][ktmp]->GetTrkFlag()>0
                         && ClusterBank[i-10][ktmp]->GetNDFlag()>0
                         && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
                         && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+10][kp],10)>0) 
                        {AlternateTriplets=true;}
                    }
                  }
                                      
                  // If everything is good, make the triplet
                  if(AlternateTriplets==false) {
                    TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-20][km],ClusterBank[i][k0],ClusterBank[i+10][kp]);
                    for(unsigned int j=0; j<seg0->GetEntries(); ++j) {seg0->GetCluster(j)->SetNDFlag(2);}

                    NDSegmentBank[i].push_back(seg0);
                  }
                }
              }
            }
          }
        }




        // Look for triplets of form m 0 X p  
        if( (i-10)>=0 && ClusterBank[i-10].size()>0 && (i+20)<=PlanesInModule && ClusterBank[i+20].size()>0 ) {

          // Look at the clusters on plane i-10...
          for(unsigned int km=0; km<ClusterBank[i-10].size(); ++km) {
              
            // ... and the clusters on plane i+20
            for(unsigned int kp=0; kp<ClusterBank[i+20].size(); ++kp) {

              if( ClusterBank[i-10][km]->GetTrkFlag()>0 && ClusterBank[i+20][kp]->GetTrkFlag()>0
                  && ClusterBank[i-10][km]->GetNDFlag()>0 && ClusterBank[i+20][kp]->GetNDFlag()>0 ) {

                // Check the track-like association of the three clusters
                TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+20][kp],10);
                   
                // If the association is good, check whether we can make the triplet
                if(TripAssocNum>0) {
                  AlternateTriplets=false;
                      
                  // Check to see if there is also an alternate triplet possibility with clusters on plane i+2
                  // i.e. Want to make sure that the X in "m 0 X p" is correct.

                  if(AlternateTriplets==false && ClusterBank[i+10].size()>0) {
                    // Look at the clusters on plane i+10
                    for(unsigned int ktmp=0; ktmp<ClusterBank[i+10].size(); ++ktmp) {
                      if(AlternateTriplets==false && ClusterBank[i+10][ktmp]->GetTrkFlag()>0 
                         && ClusterBank[i+10][ktmp]->GetNDFlag()>0
                         && ClusterBank[i+10][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0 
                         && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][km],ClusterBank[i+10][ktmp],10)>0) 
                        {AlternateTriplets=true;}
                    }
                  }
                                      
                  // If everything is good, make the triplet
                  if(AlternateTriplets==false) {
                    // Store segment on empty plane too
                    for(int j=0; j<2; ++j) {
                      TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-10][km],ClusterBank[i][k0],ClusterBank[i+20][kp]);
                      for(unsigned int k=0; k<seg0->GetEntries(); ++k) {seg0->GetCluster(k)->SetNDFlag(2);}
  
                      NDSegmentBank[i+10*j].push_back(seg0);
                    }
                  }
                }
              }
            }
          }
        }




        // Look for triplets of form m X 0 X p
        if( (i-20)>=0 && ClusterBank[i-20].size()>0 && (i+20)<=PlanesInModule && ClusterBank[i+20].size()>0 ) {

          // Look at clusters on planes i-20 and i+20
          for(unsigned int km=0; km<ClusterBank[i-20].size(); ++km) {
            for(unsigned int kp=0; kp<ClusterBank[i+20].size(); ++kp) {

              if( ClusterBank[i-20][km]->GetTrkFlag()>0 && ClusterBank[i+20][kp]->GetTrkFlag()>0
                  && ClusterBank[i-20][km]->GetNDFlag()>0 && ClusterBank[i+20][kp]->GetNDFlag()>0 ) {

                // Check the track-like association of the three clusters
                TripAssocNum=ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+20][kp],10);
                  
                // If the association is good, check whether we can make the triplet
                if(TripAssocNum>0) {
                  AlternateTriplets=false;
                    
                  // Check to see if there are alternate triplet possibilities with clusters on plane i-10 or i+10
                  // i.e. Want to make sure that the Xs in "m X 0 X p" are correct.
                    
                  // Check first X, plane i-10
                  if(AlternateTriplets==false && ClusterBank[i-10].size()>0) {
                    // Look at any clusters on plane i-10
                    for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
                      if(AlternateTriplets==false && ClusterBank[i-10][ktmp]->GetTrkFlag()>0 
                         && ClusterBank[i-10][ktmp]->GetNDFlag()>0
                         && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
                         && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+20][kp],10)>0) 
                        {AlternateTriplets=true;}
                    }
                  }
                    
                  // Check second X, plane i+10
                  if(AlternateTriplets==false && ClusterBank[i+10].size()>0) {
                    // Look at any clusters on plane i+10
                    for(unsigned int ktmp=0; ktmp<ClusterBank[i+10].size(); ++ktmp) {
                      if(AlternateTriplets==false  && ClusterBank[i+10][ktmp]->GetTrkFlag()>0 
                         && ClusterBank[i+10][ktmp]->GetNDFlag()>0
                         && ClusterBank[i+10][ktmp]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0 
                         && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i+10][ktmp],10)>0) 
                        {AlternateTriplets=true;}
                    }
                  }
                    
                  // Check both Xs together, planes i-10 and i+10
                  if(AlternateTriplets==false && ClusterBank[i-10].size()>0 
                     && ClusterBank[i+10].size()>0) {
                      
                    // Loop again at any clusters on plane i-10
                    for(unsigned int ktmp=0; ktmp<ClusterBank[i-10].size(); ++ktmp) {
                        
                      // Look again at any clusters on plane i+10
                      for(unsigned int ktmp1=0; ktmp1<ClusterBank[i+10].size(); ++ktmp1) {
                        if(AlternateTriplets==false 
                           && ClusterBank[i-10][ktmp]->GetTrkFlag()>0 && ClusterBank[i+10][ktmp1]->GetTrkFlag()>0
                           && ClusterBank[i-10][ktmp]->GetNDFlag()>0  && ClusterBank[i+10][ktmp1]->GetNDFlag()>0 
                           && ClusterBank[i-10][ktmp]->IsTrkAssoc(ClusterBank[i-20][km],ClusterBank[i][k0],10)>0 
                           && ClusterBank[i][k0]->IsTrkAssoc(ClusterBank[i-10][ktmp],ClusterBank[i+10][ktmp1],10)>0 
                           && ClusterBank[i+10][ktmp1]->IsTrkAssoc(ClusterBank[i][k0],ClusterBank[i+20][kp],10)>0)
                          {AlternateTriplets=true;}
                      }
                    }
                  }

                  // If everything is good, make the triplet
                  if(AlternateTriplets==false) {
                    // Store segment on empty planes too
                    for(int j=0; j<2; ++j) {
                      TrackSegmentCam* seg0 = new TrackSegmentCam(ClusterBank[i-20][km],ClusterBank[i][k0],ClusterBank[i+20][kp]);
                      for(unsigned int k=0; k<seg0->GetEntries(); ++k) {seg0->GetCluster(k)->SetNDFlag(2);}
        
                      NDSegmentBank[i+10*j].push_back(seg0);
                    }
                  }
                }
              }
            }
          }
        }

      }
    }
  }



  // Print out list of triplets
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ND TRIPLETS *** " << endl;
  for(int i=0; i<500; ++i) {
    for(unsigned int j=0;j<NDSegmentBank[i].size(); ++j) {
      MSG("AlgTrackCamList", Msg::kVerbose) << " Segment Num " << j << ", Plane " << i <<endl;
      
      ClusterCam* clr0 = NDSegmentBank[i][j]->GetCluster(0);
      ClusterCam* clr1 = NDSegmentBank[i][j]->GetCluster(1);
      ClusterCam* clr2 = NDSegmentBank[i][j]->GetCluster(2);
      
      MSG("AlgTrackCamList", Msg::kVerbose) << " (" << clr0->GetPlane() << "," 
                                            << clr0->GetBegTPos() << "->" << clr0->GetEndTPos() << ") "
                                            << " (" << clr1->GetPlane() << "," << clr1->GetBegTPos() 
                                            << "->" << clr1->GetEndTPos() << ") "
                                            << " (" << clr2->GetPlane() << "," << clr2->GetBegTPos() 
                                            << "->" << clr2->GetEndTPos() << ") " << endl;
    }
  }




  // Now loop over these segments and check associations
  for(int i=10; i<PlanesInModule-9; ++i) {
    
    if((i+10)<=PlanesInModule && NDSegmentBank[i].size()>0 && NDSegmentBank[i+10].size()>0) {
        
      // Look at the segments on plane i
      for(unsigned int k0=0; k0<NDSegmentBank[i].size(); ++k0) {
        TrackSegmentCam* Seg0 = NDSegmentBank[i][k0];

        // For each of these, look at the segments on plane i+10
        for(unsigned int kp=0; kp<NDSegmentBank[i+10].size(); ++kp) {
          TrackSegmentCam* NextSeg = NDSegmentBank[i+10][kp];

          // Check associations. 
          if(Seg0->IsAssoc(NextSeg)) {
              
            // If segments are associated, add to list of all associations
            Seg0->AddAssocSegToEnd(NextSeg);
            NextSeg->AddAssocSegToBeg(Seg0);
          }
        }
      }
    }
  }
  


  // Print out list of all associations
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF ALL ND TRIPLET ASSOCIATIONS *** " << endl;

  // Print out list of triplets
  for(int i=0; i<500; ++i) {
    for(unsigned int j=0;j<NDSegmentBank[i].size(); ++j) {
      TrackSegmentCam* segment = NDSegmentBank[i][j];

      MSG("AlgTrackCamList", Msg::kVerbose) << "--- Plane " << i << ", Seg Number " << j 
                                            << ", BegPlane " << segment->GetBegPlane() 
                                            << ", EndPlane " << segment->GetEndPlane() 
                                            << ", BegTPos " << segment->GetBegTPos() 
                                            << ", EndTPos " << segment->GetEndTPos() 
                                            <<endl;
      
      MSG("AlgTrackCamList", Msg::kVerbose) << " Beg: " <<endl;
      for(unsigned int k=0; k<segment->GetNAssocSegBeg(); ++k) {
        TrackSegmentCam* segtmp = segment->GetAssocSegBeg(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;
      }
      
      MSG("AlgTrackCamList", Msg::kVerbose) << " End: " <<endl;
      for(unsigned int k=0; k<segment->GetNAssocSegEnd(); ++k) {
        TrackSegmentCam* segtmp = segment->GetAssocSegEnd(k);
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Assoc number=" << k
                                              << "  begpln: " << segtmp->GetBegPlane()
                                              << "  begtpos: " << segtmp->GetBegTPos()
                                              << "  endpln: " << segtmp->GetEndPlane()
                                              << "  endtpos: " << segtmp->GetEndTPos() << endl;           
      }
    }
  }



  // Look for chains of associations and make ND 2D tracks
  bool JoinFlag;
  
  for(int i=10; i<PlanesInModule-9; ++i) {
    // Loop over the triplets centered on plane i
    for(unsigned int j=0; j<NDSegmentBank[i].size(); ++j) {

      TrackSegmentCam* Seg = NDSegmentBank[i][j];


      // Firstly, consider isolated segments. These must be their own seed segments.
      if(Seg->GetNAssocSegBeg()==0 && Seg->GetNAssocSegEnd()==0) {
        NDViewSegBank[Seg->GetPlaneView()].push_back(Seg);
        Seg->SetSeedSegment(Seg);
      }


      // Then, consider those segments which have some associations
      if(Seg->GetNAssocSegBeg()>0 || Seg->GetNAssocSegEnd()>0) {
        JoinFlag=false;

        // If there is one assoc at beginning (always false for first segment in module)
        if(Seg->GetNAssocSegBeg()==1) {
          TrackSegmentCam* Segm = Seg->GetAssocSegBeg(0);

          // Check that we are considering a simple chain of segments
          if(Segm->GetNAssocSegEnd()==1) {
            // Get the last segment in the chain of segments
            TrackSegmentCam* SeedSeg = Segm->GetSeedSegment();

            SeedSeg->AddSegment(Seg); 
            Seg->SetSeedSegment(SeedSeg);
            Seg->SetUID(-1);
            JoinFlag=true;
          }
        }

        // Now consider those triplets which aren't just one assoc at beg and one at end
        if(JoinFlag==false) {
          Seg->SetSeedSegment(Seg);
          NDViewSegBank[Seg->GetPlaneView()].push_back(Seg);
                  
          // Loop over the associations at beginning
          for(unsigned int k=0; k<Seg->GetNAssocSegBeg(); ++k) {
            TrackSegmentCam* SeedSeg = Seg->GetAssocSegBeg(k)->GetSeedSegment();
            
            Seg->AddMatchSegToBeg(SeedSeg);
            SeedSeg->AddMatchSegToEnd(Seg);
          }
        }
      }

    } 
  }




  // Print out the list of 2D tracks
  MSG("AlgTrackCamList", Msg::kVerbose)  << " *** LIST OF ND 2D TRACKS *** " << endl;
  for(int View=0; View<2; ++View) {

    MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
    for(unsigned i=0; i<NDViewSegBank[View].size(); ++i) {
      TrackSegmentCam* Seg = NDViewSegBank[View][i];
      MSG("AlgTrackCamList", Msg::kVerbose) << "  i=" << i
                                            << "  begpln: "  << Seg->GetBegPlane()
                                            << "  begtpos: " << Seg->GetBegTPos()
                                            << "  endpln: "  << Seg->GetEndPlane()
                                            << "  endtpos: " << Seg->GetEndTPos() << endl;
    }
  }




  // Now match to existing 2D tracks
  // If no match, store separately in ViewSegBank
  bool Match;

  // For ND U segments
  for(unsigned int i=0; i<NDViewSegBank[0].size(); ++i) {
    TrackSegmentCam* NDSegu = NDViewSegBank[0][i];
    Match=false;
      
    // Consider the other U segments
    for(unsigned int j=0; j<ViewSegBank[0].size(); ++j) {
      TrackSegmentCam* Segu = ViewSegBank[0][j];

      if(Segu->GetUID()>1) {

        if(NDSegu->GetEndPlane()<Segu->GetEndPlane()) {Match=NDSegu->IsAssoc(Segu);}
        else {Match=Segu->IsAssoc(NDSegu);}

        if(Match) {
          if(NDSegu->GetBegPlane()<Segu->GetBegPlane()) {Segu->AddMatchSegToBeg(NDSegu);}
          else {Segu->AddMatchSegToEnd(NDSegu);}

          MSG("AlgTrackCamList", Msg::kVerbose) << "U Match, NDbegplane " << NDSegu->GetBegPlane() 
                                                << " Endplane " << NDSegu->GetEndPlane()
                                                << " Begtpos " << NDSegu->GetBegTPos() 
                                                << " Endtpos " << NDSegu->GetEndTPos()
                                                << " 2d track begplane " << Segu->GetBegPlane() 
                                                << " Endplane " << Segu->GetEndPlane()
                                                << " Begtpos " << Segu->GetBegTPos() 
                                                << " Endtpos " << Segu->GetEndTPos() << endl;
        }
      }
    }

    if(Match) {NDSegu->SetUID(2);}
    else {NDSegu->SetUID(2); TempNDViewSegBank[0].push_back(NDSegu);}
  }


  // For ND V segments
  for(unsigned int i=0; i<NDViewSegBank[1].size(); ++i) {
    TrackSegmentCam* NDSegv = NDViewSegBank[1][i];
    Match=false;
        
    // Consider the other V segments
    for(unsigned int j=0; j<ViewSegBank[1].size(); ++j) {
      TrackSegmentCam* Segv = ViewSegBank[1][j];

      if(Segv->GetUID()>1) {

        if(NDSegv->GetEndPlane()<Segv->GetEndPlane()) {Match=NDSegv->IsAssoc(Segv);}
        else {Match=Segv->IsAssoc(NDSegv);}

        if(Match) {
          if(NDSegv->GetBegPlane()<Segv->GetBegPlane()) {Segv->AddMatchSegToBeg(NDSegv);}
          else {Segv->AddMatchSegToEnd(NDSegv);}

          MSG("AlgTrackCamList", Msg::kVerbose) << "V Match, NDbegplane " << NDSegv->GetBegPlane() 
                                                << " Endplane " << NDSegv->GetEndPlane()
                                                << " Begtpos " << NDSegv->GetBegTPos() 
                                                << " Endtpos " << NDSegv->GetEndTPos()
                                                << " 2d track begplane " << Segv->GetBegPlane() 
                                                << " Endplane " << Segv->GetEndPlane()
                                                << " Begtpos " << Segv->GetBegTPos() 
                                                << " Endtpos " << Segv->GetEndTPos() << endl;
        }
      }
    }

    if(Match) {NDSegv->SetUID(2);}
    else {NDSegv->SetUID(2); TempNDViewSegBank[1].push_back(NDSegv);}
  }


  // Now store those segments marked above in ViewSegBank
  for(int View=0; View<2; ++View) {
    for(unsigned int i=0; i<TempNDViewSegBank[View].size(); ++i) {
      ViewSegBank[View].push_back(TempNDViewSegBank[View][i]);
    }
  }


  // Clean up
  for(int i=0; i<2; ++i) {NDViewSegBank[i].clear(); TempNDViewSegBank[i].clear();}

  return;
}



void AlgTrackCamList::Form2DTracks()
{
  // Of the segments identified as good above, we identify the segments which
  // are good seeds for the track i.e. from which we can propagate back and forth
  // along matched segments to find a long track.
 
  // Then, for the first seed segment, we try to propagate backwards and forwards,
  // marking the segments we use with different TmpTrkFlag settings.

  // For a seed segment, Seg0, we firstly move in the direction of increasing plane
  // number. We mark all the segments we use with TmpTrkFlag 1. The seed segment is
  // labelled with TmpTrkFlag 3.

  // Some paths will lead further than others. We are interested in the longest paths, 
  // and move back towards Seg0 along the longest paths, marking the segments used 
  // with TmpTrkFlag 2.

  // Having done this, we can carry out a similar procedure, but initially moving 
  // backwards from Seg0 to lower plane numbers, before returning to Seg0 again.

  // So, this arrangement of segments...
  //
  //       Seg  Seg                                  Seg
  //                  Seg                      Seg
  //                        Seg   Seg0   Seg 
  //                  Seg                      Seg   Seg   Seg
  //                                                 
  //               Seg                               Seg
  //                                                       Seg   Seg

  // Is labelled as follows...
  //
  //       2     2                                    1 
  //                   2                        1 
  //                         2     3      2  
  //                   1                        2     2     2 
  //                                                 
  //                1                                 2 
  //                                                        2     2 
  //

  // It is possible that there are multiple paths along segments with TmpTrkFlag 2, 
  // as above. The best path is selected by obtaining a score for each, based on
  // straightness and length.



  // Define containers
  vector<TrackSegmentCam*> BestSeedSegments;
  vector<TrackSegmentCam*> SeedSegments;

  vector<TrackSegmentCam*> Temp;
  vector<TrackSegmentCam*> mTemp;
  vector<TrackSegmentCam*> pTemp;
  vector<TrackSegmentCam*> BegTemp;
  vector<TrackSegmentCam*> EndTemp;
  vector<TrackSegmentCam*> TempBank1;
  vector<TrackSegmentCam*> TempBank2;

  vector <vector<TrackSegmentCam*> > BegBank;
  vector <vector<TrackSegmentCam*> > EndBank;

  vector <vector<TrackSegmentCam*> > TempBegBank;
  vector <vector<TrackSegmentCam*> > TempEndBank;


  bool Cont; bool tmpFlag; bool TrackFlag;
  int ntrks; int Counter; int nplane; int npts;
  int ObjCounter; int ObjVectorCounter; const int fObjVectorMax=100000;
  int id; double Score; double TopScore;

  unsigned int MostClusters;
  bool AlreadyAdded;


  // Main loop is over the views, first for U, then V
  if(ViewSegBank[0].size()>0 && ViewSegBank[1].size()>0) {
    for(int View=0; View<2; ++View) {
    
      Cont=true; ntrks=0;

      while(Cont==true) {
        Cont=false;
        BegBank.clear(); EndBank.clear();



        // First selection of SEED SEGMENTS - segment from which we can move
        // along matched associations to find majority of the track.
        //
        // First guesses at seed segments are stored in SeedSegments Container
        BestSeedSegments.clear(); SeedSegments.clear();

        // Loop over the entries in ViewSegBank
        for(unsigned int i=0; i<ViewSegBank[View].size(); ++i) {
          TrackSegmentCam* Seg = ViewSegBank[View][i];

          if(Seg->GetUID()>0 && Seg->GetUID()<3) {
            Counter=0;
            
            // Loop over the constituent clusters and count the number that are track-like
            for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
              ClusterCam* Clust = Seg->GetCluster(j);
              if(Clust->GetTrkPlnFlag()>0 && Clust->GetTrkFlag()==2) {Counter++;}
            }
            
            // If we've already made a track and this segment doesn't contain 3 track-like
            // clusters, don't consider it any further
            if(ntrks>0 && Counter<3) {Seg->SetUID(0);}
          }

          // Temporarily store all the likely seed segments. Those with UID 2 are most track-like.
          if(Seg->GetUID()==2) {BestSeedSegments.push_back(Seg);}
          if(Seg->GetUID()==1) {SeedSegments.push_back(Seg);}
        }


        // Hopefully we have segments with UID 2 (most track-like), else take those with UID 1
        // (just track-like) and move them from SeedSegments to BestSeedSegments. 
        // Can then empty SeedSegments.
        if(BestSeedSegments.size()==0) {BestSeedSegments=SeedSegments;}
        
        SeedSegments.clear();

        SeedSegments=BestSeedSegments;




        // Now, using our initial seed segments, try propagating back and forth
        // from each segment to refine our choice.

        // Final SEED SEGMENTS are stored BestSeedSegments Container
        BestSeedSegments.clear();


        // If we have some seed segments to work with, loop over them
        if(SeedSegments.size()>0) { 
          nplane=-1;

          for(unsigned int i=0; i<SeedSegments.size(); ++i) {

            // See how far we can propagate backwards to lower plane numbers
            mTemp.clear(); BegTemp.clear(); 

            // First put seed segment into BegTemp then begin while loop
            BegTemp.push_back(SeedSegments[i]); TempBank1.clear(); 

            while(BegTemp.size()>0) {
              Temp.clear();
              
              // Copy any segments into a temporary holder and loop over them,
              // gradually propagating back to lower plane numbers
              Temp=BegTemp; BegTemp.clear();
            
              for(unsigned int j=0; j<Temp.size(); ++j) {
                tmpFlag=false;
                
                // For each segment, loop over their beginning matches
                for(unsigned int k=0; k<Temp[j]->GetNMatchSegBeg(); ++k) {
                  TrackSegmentCam* Segb = Temp[j]->GetMatchSegBeg(k);
                  
                  // If matched segment is track-like and has zero TmpTrkFlag, 
                  // push it back into BegTemp, so the while loop continues.
                  if(Segb->GetUID()>-1 && Segb->GetUID()<3) {
                    tmpFlag=true;
                    
                    if(Segb->GetTmpTrkFlag()<1) {
                      BegTemp.push_back(Segb);
                      Segb->SetTmpTrkFlag(1); // Temporarily mark segments used
                      TempBank1.push_back(Segb);
                    }
                  }
                }
                
                // Must have gone as far as we can go with this seed segment. Store segment at
                // lowest plane.
                if(tmpFlag==false) {mTemp.push_back(Temp[j]);}
              }
            } // End while BegTemp.size()>0

            // Set the flags back to zero, for when we carry out actual propagation, below.
            for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}


            
            // Now see how far we can propagate forwards to higher plane numbers
            pTemp.clear(); EndTemp.clear();

            // Put the seed segment into EndTemp then start while loop
            EndTemp.push_back(SeedSegments[i]); TempBank1.clear();
            
            while(EndTemp.size()>0) {
              Temp.clear();
              
              // Copy any segments into a temporary holder and loop over them,
              // gradually propagating to higher plane numbers
              Temp=EndTemp; EndTemp.clear();
             
              for(unsigned int j=0; j<Temp.size(); ++j) {
                tmpFlag=false;
                
                // For each segment, loop over their end matches
                for(unsigned int k=0; k<Temp[j]->GetNMatchSegEnd(); ++k) {
                  TrackSegmentCam* Sege = Temp[j]->GetMatchSegEnd(k);

                  // If the matched segment is track-like and has zero TmpTrkFlag, 
                  // push back into EndTemp, so the while loop continues.
                  if(Sege->GetUID()>-1 && Sege->GetUID()<3) {
                    tmpFlag=true;
                    
                    if(Sege->GetTmpTrkFlag()<1) {
                      EndTemp.push_back(Sege);
                      Sege->SetTmpTrkFlag(1); // Temporarily mark segments used
                      TempBank1.push_back(Sege);
                    }
                  }
                }
                
                // Must have gone as far as we can go with this seed segment. Store segment at
                // highest plane.
                if(tmpFlag==false) {pTemp.push_back(Temp[j]);}
              }
            } // End while EndTemp.size()>0
            
            // Set the flags back to zero, for when we carry out actual propagation, below.
            for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}
            
            
            
            // Find the maximum number of planes we can span by propagating the seed segment
            // back and forth.
            npts=-1;
            for(unsigned int j=0; j<mTemp.size(); ++j) {
              for(unsigned int k=0; k<pTemp.size(); ++k) {
                if((1+pTemp[k]->GetEndPlane()-mTemp[j]->GetBegPlane())>npts) {
                  npts=1+pTemp[k]->GetEndPlane()-mTemp[j]->GetBegPlane();
                }
              }
            }
            
            SeedSegments[i]->SetNPlanes(npts);
            if(npts>nplane) {nplane=npts;}
          
          } // End loop over seed segments stored above


          // Store the seed segments that lead to the biggest propagation in BestSeedSegments
          if(nplane>0) {
            for(unsigned int i=0; i<SeedSegments.size(); ++i) {
              if(SeedSegments[i]->GetNPlanes()>nplane-5) {BestSeedSegments.push_back(SeedSegments[i]);}       
            }
          }
          
          
        } // End finding final seed segments
        

        // Order contents of BestSeedSegments
        SeedSegments=BestSeedSegments; BestSeedSegments.clear();
        
        for(unsigned int i=0; i<SeedSegments.size(); ++i) {
          MostClusters=0;
          TrackSegmentCam* LargestSeg = 0;        
          
          for(unsigned int j=0; j<SeedSegments.size(); ++j) {
            AlreadyAdded=false;
            
            for(unsigned int k=0; k<BestSeedSegments.size(); ++k) {
              if(BestSeedSegments[k]==SeedSegments[j]) {AlreadyAdded=true;}
            }
            if(AlreadyAdded) {continue;}

            if(SeedSegments[j]->GetEntries()>MostClusters) {
              MostClusters=SeedSegments[j]->GetEntries();
              LargestSeg=SeedSegments[j];
            }
          }

          if(LargestSeg) {BestSeedSegments.push_back(LargestSeg);}
        }

        

        // Now use the final seed segments to actually PROPAGATE back and forth, 
        // setting the TmpTrkFlags for the segments used in the propagation.
        TempBank2.clear();

        // Loop over the seed segments, which are stored in BestSeedSegments
        for(unsigned int i=0; i<BestSeedSegments.size(); ++i) {
          TrackSegmentCam* Seg0 = BestSeedSegments[i];

          MSG("AlgTrackCamList", Msg::kVerbose) << " Seed BegPlane" << Seg0->GetBegPlane() 
                                                << " Endplane " << Seg0->GetEndPlane()
                                                << " Begtpos " << Seg0->GetBegTPos() 
                                                << " Endtpos " << Seg0->GetEndTPos()
                                                << " UID " << Seg0->GetUID()
                                                << " entries " << Seg0->GetEntries() << endl;


          // If we make a track, we will definitely include the seed segment, so set
          // it's TmpTrkFlag to 3.
          Seg0->SetTmpTrkFlag(3); 
          TrackFlag=false;

          TempBank1.clear(); TempBank1.push_back(Seg0); 


          // Track backwards (1)
          mTemp.clear(); BegTemp.clear();

          // Push the seed segment into BegTemp and begin while loop
          BegTemp.push_back(Seg0);

          while(BegTemp.size()>0) {
            Temp.clear();

            // Copy any segments into a temporary holder and loop over them,
            // gradually propagating back to lower plane numbers
            Temp=BegTemp; BegTemp.clear();

            for(unsigned int j=0; j<Temp.size(); ++j) {
              tmpFlag=false;
              TrackSegmentCam* Segtmp = Temp[j];
              
              // For each segment, loop over their beginning matches
              for(unsigned int k=0; k<Segtmp->GetNMatchSegBeg(); ++k){
                TrackSegmentCam* Segb = Segtmp->GetMatchSegBeg(k);
                
                // If the matched segment is track-like and has zero TmpTrkFlag, 
                // push back into BegTemp, so the while loop continues.
                if(Segb->GetUID()>-1 && Segb->GetUID()<3) {
                  tmpFlag=true;
                  
                  if(Segb->GetTmpTrkFlag()<1) {
                    BegTemp.push_back(Segb);
                    Segb->SetTmpTrkFlag(1); // Mark the segments with TmpTrkFlag 1
                    TempBank1.push_back(Segb);
                  }
                }
              }

              // Must have gone as far as we can go with this seed segment. Store segment at
              // lowest plane.
              if(tmpFlag==false) {mTemp.push_back(Segtmp);}
            }
          }
          
          // Note the first plane of the potential 2D track
          nplane=999;
          for(unsigned int j=0; j<mTemp.size(); ++j) {
            if(mTemp[j]->GetBegPlane()<nplane) {nplane=mTemp[j]->GetBegPlane();}
          }


          // For the longest possible paths, move back towards the seed segment,
          // setting TmpTrkFlags to 2 for the segments used.

          // Find the segments reached by the longest propagations
          BegTemp.clear();
          for(unsigned int j=0; j<mTemp.size(); ++j) {
            if(mTemp[j]->GetBegPlane()<nplane+5) {BegTemp.push_back(mTemp[j]);}
          }


          // From these segments, propagate forwards and set the flags.
          while(BegTemp.size()>0) {
            Temp.clear();
            
            for(unsigned int j=0; j<BegTemp.size(); ++j) {
              if(BegTemp[j]->GetTmpTrkFlag()<2) {
                BegTemp[j]->SetTmpTrkFlag(2); // Mark the segments with TmpTrkFlag 2
                Temp.push_back(BegTemp[j]);
              }
              
              if(BegTemp[j]->GetTrkFlag()<1) {
                BegTemp[j]->SetTrkFlag(1);    // Also mark the segments with TrkFlag 1
                TempBank2.push_back(BegTemp[j]);
                TrackFlag=true;
              }
            }
            BegTemp.clear();
            
            
            for(unsigned int j=0; j<Temp.size(); ++j) {
              for(unsigned int k=0; k<Temp[j]->GetNMatchSegEnd(); ++k) {
                if(Temp[j]->GetMatchSegEnd(k)->GetTmpTrkFlag()>0) {
                  BegTemp.push_back(Temp[j]->GetMatchSegEnd(k));
                }
              }
            }
          }
          


          // Track forwards (1)
          pTemp.clear(); EndTemp.clear();
          
          
          // Push seed segment into EndTemp and begin while loop
          EndTemp.push_back(Seg0);

          while(EndTemp.size()>0) {
            Temp.clear();
         
            // Copy any segments into a temporary holder and loop over them,
            // gradually propagating to higher plane numbers
            Temp=EndTemp; EndTemp.clear();

            for(unsigned int j=0; j<Temp.size(); ++j) {
              tmpFlag=false;
              TrackSegmentCam* Segtmp = Temp[j];
              
              // For each segment, loop over their end matches
              for(unsigned int k=0; k<Segtmp->GetNMatchSegEnd(); ++k){
                TrackSegmentCam* Sege = Segtmp->GetMatchSegEnd(k);
                
                // If the matched segment is track-like and has zero TmpTrkFlag, 
                // push back into BegTemp, so the while loop continues.
                if(Sege->GetUID()>-1 && Sege->GetUID()<3) {
                  tmpFlag=true;
                  
                  if(Sege->GetTmpTrkFlag()<1) {
                    EndTemp.push_back(Sege);
                    Sege->SetTmpTrkFlag(1); // Mark the segments with TmpTrkFlag 1
                    TempBank1.push_back(Sege);
                  }
                }
              }

              // Must have gone as far as we can go with this seed segment. Store segment at
              // highest plane.
              if(tmpFlag==false) {pTemp.push_back(Segtmp);}
            }
          }

          // Note the end plane of potential 2D track     
          nplane=-999;
          for(unsigned int j=0; j<pTemp.size(); ++j) {
            if(pTemp[j]->GetEndPlane()>nplane) {nplane=pTemp[j]->GetEndPlane();}
          }


          // For the longest possible paths, move back towards the seed segment,
          // setting TmpTrkFlags to 2 for the segments used.

          // Find the segments reached by the longest propagations
          EndTemp.clear();
          for(unsigned int j=0; j<pTemp.size(); ++j) {
            if(pTemp[j]->GetEndPlane()>nplane-5) {EndTemp.push_back(pTemp[j]);}
          }

          
          // From these segments, propagate back and set the flags.
          while(EndTemp.size()>0) {
            Temp.clear();
            
            for(unsigned int j=0; j<EndTemp.size(); ++j) {
              if(EndTemp[j]->GetTmpTrkFlag()<2) {
                EndTemp[j]->SetTmpTrkFlag(2); // Mark the segments with TmpTrkFlag 2
                Temp.push_back(EndTemp[j]);
              }
              
              if(EndTemp[j]->GetTrkFlag()<1) {
                EndTemp[j]->SetTrkFlag(1);    // Also mark the segments with TrkFlag 1
                TempBank2.push_back(EndTemp[j]);
                TrackFlag=true;
              }
            }
            EndTemp.clear();
            
            for(unsigned int j=0; j<Temp.size(); ++j) {
              for(unsigned int k=0; k<Temp[j]->GetNMatchSegBeg(); ++k) {
                if(Temp[j]->GetMatchSegBeg(k)->GetTmpTrkFlag()>0) {
                  EndTemp.push_back(Temp[j]->GetMatchSegBeg(k));
                }
              }
            }
          }

          // End setting TmpTrkFlags




          // From paths of segments with TmpTrkFlag==2, find the best set of
          // segments at the beginning and the best set of segments at the end

          // If propagation from the seed segment has given us a good potential 2D track
          MSG("AlgTrackCamList", Msg::kVerbose) << "TrackFlag " << TrackFlag << endl;

          if(TrackFlag==true) { 

            // Track backwards (2) and find best set of beginning segments
            TempBegBank.clear(); ObjVectorCounter=0;
            vector<TrackSegmentCam*> TmpBegSegBank;

            TmpBegSegBank.push_back(Seg0);
            TempBegBank.push_back(TmpBegSegBank);
            ObjCounter=1;
            

            // Store the series of beginning segments that were flagged above during the propagation.
            // Each entry in TempBegbank will be a vectors of segments (a series of beginning segments 
            // flagged above).

            // TempBegBank will contain vectors of segments representing every possible path back through
            // the segments marked with TmpTrkFlag 2 (possible beginning sections for the 2D track)
            while(ObjCounter>0 && ObjVectorCounter<fObjVectorMax) {
              ObjCounter=0;
              npts=TempBegBank.size();

              for(int j=0; j<npts; ++j) {
                TrackSegmentCam* Segtmp = TempBegBank[j].back();
                mTemp.clear();
                
                for(unsigned int k=0; k<Segtmp->GetNMatchSegBeg(); ++k) {
                  TrackSegmentCam* Segbeg = Segtmp->GetMatchSegBeg(k);
                  if(Segbeg->GetTmpTrkFlag()>1) {mTemp.push_back(Segbeg);};
                }

                if(mTemp.size()>0) {
                  for(unsigned int k=1; k<mTemp.size(); ++k) {
                    if(ObjVectorCounter<fObjVectorMax) { 
                      vector<TrackSegmentCam*> NewObj = TempBegBank[j];
                      
                      NewObj.push_back(mTemp[k]);
                      TempBegBank.push_back(NewObj);

                      ObjVectorCounter++;
                    }
                    else {MSG("AlgTrackCamList", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl; break;}
                  }

                  if(ObjVectorCounter>=fObjVectorMax) {break;}

                  TempBegBank[j].push_back(mTemp[0]); ObjCounter++;
                }
              }
            } // End while ObjCounter>0

            // Find beginning section with best score and store it
            MSG("AlgTrackCamList", Msg::kVerbose) << "beginning score " << endl;

            id=-1; TopScore=-1.;
            for(unsigned int j=0; j<TempBegBank.size(); ++j) {
              Score = Seg0->GetScore(&TempBegBank[j],0);
              if(Score>TopScore) {TopScore=Score; id=j;}
            }

            if(id!=-1) {BegBank.push_back(TempBegBank[id]);}
            TempBegBank.clear();




            // Track forwards (2) and find best set of end segments
            TempEndBank.clear(); ObjVectorCounter=0;
            vector<TrackSegmentCam*> TmpEndSegBank;

            TmpEndSegBank.push_back(Seg0);
            TempEndBank.push_back(TmpEndSegBank);
            ObjCounter=1;
            

            // Store the series of end segments that were flagged above during the propagation.
            // Each entry in TempEndbank will be a vectors of segments (a series of end segments 
            // flagged above).

            // TempBegBank will contain vectors of segments representing every possible path forward through
            // the segments marked with TmpTrkFlag 2 (possible end sections for the 2D track)
            while(ObjCounter>0 && ObjVectorCounter<fObjVectorMax) {
              ObjCounter=0;
              npts=TempEndBank.size();
              
              for(int j=0; j<npts; ++j) {
                TrackSegmentCam* Segtmp = TempEndBank[j].back();

                pTemp.clear();
                
                for(unsigned int k=0; k<Segtmp->GetNMatchSegEnd(); ++k) {
                  TrackSegmentCam* Segend = Segtmp->GetMatchSegEnd(k);
                  if(Segend->GetTmpTrkFlag()>1) {pTemp.push_back(Segend);};
                }

                if(pTemp.size()>0) {
                  for(unsigned int k=1; k<pTemp.size(); ++k) {
                    if(ObjVectorCounter<fObjVectorMax) {
                      vector<TrackSegmentCam*> NewObj = TempEndBank[j];
                      
                      NewObj.push_back(pTemp[k]);
                      TempEndBank.push_back(NewObj);

                      ObjVectorCounter++;
                    }
                    else {MSG("AlgTrackCamList", Msg::kWarning) << " *** CRAZY MEMORY GUZZLING *** " << endl; break;}
                  }

                  if(ObjVectorCounter>=fObjVectorMax) {break;}

                  TempEndBank[j].push_back(pTemp[0]); ObjCounter++;
                }
              }
            } // End while ObjCounter>0

            // Find end section with best score and store it
            MSG("AlgTrackCamList", Msg::kVerbose) << "end score " << endl;

            id=-1; TopScore=-1.;
            for(unsigned int j=0; j<TempEndBank.size(); ++j) {
              Score = Seg0->GetScore(0,&TempEndBank[j]);
              if(Score>TopScore) {TopScore=Score; id=j;}
            }

            if(id!=-1) {EndBank.push_back(TempEndBank[id]);}
            TempEndBank.clear(); // maybe clear individual vectors of segments
          }

          
          // Clear up
          for(unsigned int j=0; j<TempBank1.size(); ++j) {TempBank1[j]->SetTmpTrkFlag(0);}

        } // End loop over BestSeedSegments
        for(unsigned int j=0; j<TempBank2.size(); ++j) {TempBank2[j]->SetTrkFlag(0);}





        // Using the best beginning sections and the best end sections from all 
        // seed segments, find the BEST COMPLETE 2D track
        if(BegBank.size()>0 && EndBank.size()>0) {

          // Calculate a score, including the seed segment and best beginning/end paths,
          // to find the 'best' complete track
          MSG("AlgTrackCamList", Msg::kVerbose) << "overall score " << endl;

          id=-1; TopScore=-1.;
          for(unsigned int i=0; i<BegBank.size(); ++i) {
            Score = BegBank[i][0]->GetScore(&BegBank[i],&EndBank[i]);
            if(Score>TopScore) {TopScore=Score; id=i;}
          }

          // If we have found a best track, add together the segments to form
          // a segment that is the 2D track, giving this UID 3. Mark segments  
          // used by setting their UID to -1.
          if(1+id>0) {
            // Get the seed segment
            TrackSegmentCam* Seg0 = BegBank[id][0];

            // Set segment UIDs to -1 if they are part of the track, so they
            // can't be used in another track.
            for(unsigned int i=1; i<BegBank[id].size(); ++i) {
              TrackSegmentCam* Segbeg = BegBank[id][i];
              Seg0->AddSegment(Segbeg);
              Segbeg->SetUID(-1);
            }

            for(unsigned int i=1; i<EndBank[id].size(); ++i) {
              TrackSegmentCam* Segend = EndBank[id][i];
              Seg0->AddSegment(Segend);
              Segend->SetUID(-1);
            }
            
            // Mark long segment as a 2D track by setting its UID to 3
            Seg0->SetUID(3);

            // Mark all the clusters in 2D track with TrkFlag 3
            for(unsigned int i=0; i<Seg0->GetEntries(); ++i) {
              Seg0->GetCluster(i)->SetTrkFlag(3);
            }

            // Store for use in Join2D tracks method
            PossibleJoins[View].push_back(Seg0);

            // Can now look for more 2D tracks
            Cont=true; ntrks++;
          }
        }

        // Clean up
        BegBank.clear(); EndBank.clear();

        
      } // End while Cont==true loop
    } // End loop over planeviews
  } // End check to see if there is anything in ViewSegBank
  

  
  // Print out the list of 2D tracks
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF 2D TRACKS *** " << endl;
  for(int View=0; View<2; ++View) {

    MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
    for(unsigned i=0; i<ViewSegBank[View].size(); ++i) {
      TrackSegmentCam* Seg = ViewSegBank[View][i];
      if(Seg->GetUID()>2) {
        MSG("AlgTrackCamList", Msg::kVerbose) << "  i=" << i
                                              << "  begpln: "  << Seg->GetBegPlane()
                                              << "  begtpos: " << Seg->GetBegTPos()
                                              << "  endpln: "  << Seg->GetEndPlane()
                                              << "  endtpos: " << Seg->GetEndTPos() << endl;
      }
    }
  }

  return;
}




void AlgTrackCamList::Join2DTracks()
{
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** JOIN 2D TRACKS *** " << endl;

  double ShortestDist; double LongestAllowedDist=5.;
  double DeltaTPos, DeltaZPos, Dist;
  int DeltaPlane, Length, UID, NumAgree;

  TrackSegmentCam* Seg1ToAdd=0;
  TrackSegmentCam* Seg2ToAdd=0;
  bool PossibleAssocs, Cont, Overlap, DiffGrad;

  if(ModuleType==2) {LongestAllowedDist=2.;}


  // Deal with two nearby tracks in ND  
  if(ModuleType==2 && PossibleJoins[0].size()==2 && PossibleJoins[1].size()==2) {
    Overlap=false;

    // Check first to see if the 2D tracks overlap and share clusters
    for(int View=0; View<2; ++View) {
      NumAgree=0; 

      for(unsigned int i=0; i<PossibleJoins[View][0]->GetEntries(); ++i) {
        ClusterCam* Clust = PossibleJoins[View][0]->GetCluster(i);
        
        for(unsigned int j=0; j<PossibleJoins[View][1]->GetEntries(); ++j) {
          if(PossibleJoins[View][1]->GetCluster(j)==Clust) {Overlap=true; break;}
        }
        if(Overlap) {break;}

      }
      if(Overlap) {break;}
    }

    

    // Check for very different gradients
    DiffGrad=false;

    for(int View=0; View<2; ++View) {
      if(PossibleJoins[View][0]->GetBegPlane()<PossibleJoins[View][1]->GetBegPlane()) {
        if(fabs(PossibleJoins[View][0]->GetEndDir()
                -PossibleJoins[View][1]->GetBegDir())>0.7) {DiffGrad=true;}
      }
      else {
        if(fabs(PossibleJoins[View][1]->GetEndDir()
                -PossibleJoins[View][0]->GetBegDir())>0.7) {DiffGrad=true;}  
      }
      if(DiffGrad) {break;}
    }

    if(DiffGrad) {
      PossibleJoins[0].clear(); 
      PossibleJoins[1].clear();
      return;
    }

    

    if(Overlap==false) {
      // Check associations apply in both views
      if( (PossibleJoins[0][0]->IsAssoc(PossibleJoins[0][1])
           || PossibleJoins[0][1]->IsAssoc(PossibleJoins[0][0]))
          && (PossibleJoins[1][0]->IsAssoc(PossibleJoins[1][1])
              || PossibleJoins[1][1]->IsAssoc(PossibleJoins[1][0])) ) {
      
        PossibleJoins[0][0]->AddSegment(PossibleJoins[0][1]); 
        PossibleJoins[0][1]->SetUID(-1);

        PossibleJoins[1][0]->AddSegment(PossibleJoins[1][1]); 
        PossibleJoins[1][1]->SetUID(-1);        
      }

      // Clear up and leave Join2DTracks method
      PossibleJoins[0].clear(); PossibleJoins[1].clear();
    
      return;
    }
  }



  else if(ModuleType==2 && PossibleJoins[0].size()>2 && PossibleJoins[1].size()>2 
          && PossibleJoins[0].size()==PossibleJoins[1].size()) {
    
    const unsigned int Size = PossibleJoins[0].size();

    bool UView[100], VView[100]; bool AllMatch=true;
    int BegPlane, EndPlane;

    if(Size<100) {
      for(unsigned int i=0; i<Size; ++i) {UView[i]=false; VView[i]=false;}

      for(unsigned int i=0; i<Size; ++i) {
        BegPlane=PossibleJoins[0][i]->GetBegPlane();
        EndPlane=PossibleJoins[0][i]->GetEndPlane();
        
        for(unsigned int j=0; j<Size; ++j) {
          
          if(VView[j]) {continue;}

          if(abs(BegPlane-PossibleJoins[1][j]->GetBegPlane())<8 
             && abs(EndPlane-PossibleJoins[1][j]->GetEndPlane())<8) {
            
            UView[i]=true; VView[j]=true; break;
          }
        }
      }
      
      for(unsigned int i=0; i<Size; ++i) {if(UView[i]==false || VView[i]==false) {AllMatch=false;}}
      
      if(AllMatch) {
        PossibleJoins[0].clear(); 
        PossibleJoins[1].clear(); 
        return;
      }
    }
  }



  // First, make the best possible joins between the tracks
  for(int View=0; View<2; ++View) {
    PossibleAssocs=true;

    while(PossibleAssocs) {
      PossibleAssocs=false; 
      ShortestDist=LongestAllowedDist; Seg1ToAdd=0; Seg2ToAdd=0; 

      // Consider all pairs of 2D tracks and find the single best association.
      for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
        TrackSegmentCam* Seg1 = PossibleJoins[View][i];

        if(Seg1->GetUID()>2) {

          for(unsigned int j=0; j<PossibleJoins[View].size(); ++j) {
            TrackSegmentCam* Seg2 = PossibleJoins[View][j];
          
            if(Seg2->GetUID()>2 && Seg1!=Seg2) {         

              // Calculate distance between segments
              DeltaTPos=fabs(Seg2->GetBegTPos()-Seg1->GetEndTPos());
              DeltaZPos=fabs(Seg2->GetBegZPos()-Seg1->GetEndZPos());
              DeltaPlane=Seg2->GetBegPlane()-Seg1->GetEndPlane();

              Dist=pow(DeltaTPos*DeltaTPos + DeltaZPos*DeltaZPos,0.5);

              if(DeltaPlane>=0 && Dist<ShortestDist)  {
                if(Seg1->IsAssoc(Seg2)) {Seg1ToAdd=Seg1; Seg2ToAdd=Seg2; ShortestDist=Dist;}
              }
              else if(DeltaPlane<=0 && Dist<ShortestDist) {
                if(Seg2->IsAssoc(Seg1)) {Seg1ToAdd=Seg1; Seg2ToAdd=Seg2; ShortestDist=Dist;}
              }
            }
          }
        }
      }


      // If a good association is found, add the segments and look for the next association.
      if(ShortestDist<LongestAllowedDist && Seg2ToAdd!=0 && Seg1ToAdd!=0) { 
        MSG("AlgTrackCamList", Msg::kVerbose) << "Missing 2D track assocation: Seg1 " 
                                              << Seg1ToAdd->GetBegPlane() << " " << Seg1ToAdd->GetEndPlane() 
                                              << " (" << Seg1ToAdd->GetBegTPos() << "," << Seg1ToAdd->GetEndTPos() 
                                              << "), Seg2 " << Seg2ToAdd->GetBegPlane() 
                                              << " " << Seg2ToAdd->GetEndPlane() << " (" << Seg2ToAdd->GetBegTPos() 
                                              << "," << Seg2ToAdd->GetEndTPos() << endl;

        Seg1ToAdd->AddSegment(Seg2ToAdd); Seg2ToAdd->SetUID(-1); PossibleAssocs=true;
      }
    }



    // Now find the longest possible 2D track and remove 2d tracks that overlap.
    Cont=true; 

    while(Cont) {
      Cont=false;

      TrackSegmentCam* BestTrk = 0;
      Length=0;

      for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
        TrackSegmentCam* Seg = PossibleJoins[View][i];
        if(Seg->GetUID()>0 && Seg->GetUID()<4 && (Seg->GetEndPlane()-Seg->GetBegPlane())>Length) {
          BestTrk=Seg; Length=Seg->GetEndPlane()-Seg->GetBegPlane();
        }
      }


      if(BestTrk) {
        BestTrk->SetUID(4);

        for(unsigned int i=0; i<BestTrk->GetEntries(); ++i) {
          ClusterCam* BestTrkClust = BestTrk->GetCluster(i);
        
          if(BestTrkClust->GetPlane()==BestTrk->GetBegPlane() || BestTrkClust->GetPlane()==BestTrk->GetEndPlane()) {continue;}

          for(unsigned int j=0; j<PossibleJoins[View].size(); ++j) {
            TrackSegmentCam* Seg = PossibleJoins[View][j];
            NumAgree=0; 

            if(Seg->GetUID()<0 || Seg->GetUID()>3) {continue;}
  
            for(unsigned int k=0; k<Seg->GetEntries(); ++k) {
              ClusterCam* SegClust = Seg->GetCluster(k);

              if(SegClust==BestTrkClust) {NumAgree++;}

              if(NumAgree>2) {Seg->SetUID(-1); break;}
            }
          }
        }
      }

      // Conditions for continuing here
      for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
        UID=PossibleJoins[View][i]->GetUID();

        if(UID>0 && UID<4) {Cont=true;}
      }
      
    }

   
    // Reset UIDs
    for(unsigned int i=0; i<PossibleJoins[View].size(); ++i) {
      if(PossibleJoins[View][i]->GetUID()==4) {PossibleJoins[View][i]->SetUID(3);}
    }

    PossibleJoins[View].clear();
  }



  // Print out the list of 2D tracks
  for(int View=0; View<2; ++View) {
    MSG("AlgTrackCamList", Msg::kVerbose) << "View: " << View << endl;
    for(unsigned i=0; i<ViewSegBank[View].size(); ++i) {
      TrackSegmentCam* Seg = ViewSegBank[View][i];
      if(Seg->GetUID()>2) {
        MSG("AlgTrackCamList", Msg::kVerbose)  << "  i=" << i
                                               << "  begpln: "  << Seg->GetBegPlane()
                                               << "  begtpos: " << Seg->GetBegTPos()
                                               << "  endpln: "  << Seg->GetEndPlane()
                                               << "  endtpos: " << Seg->GetEndTPos() << endl;
      }
    }
  }

  return;
}




void AlgTrackCamList::SecondUVComparison()
{
  // Look at the 2D tracks formed and try to match 2D tracks in the U view
  // with those in the V view, based on a rough overlap in plane numbers and time.
  // Store those that match.

  // 2D tracks are only ever formed within a single SuperModule. At this point
  // we try to match FD 2D tracks across the SM gap. Using the 2D tracks we have
  // just stored, we check association across the gap. Any final ambiguity is
  // resolved by getting a score for each possible match. 

  // If things look good, U/V partners are set and 2D tracks are joined across the gap.

  // Any other possible 2D track matches between views are then considered and the
  // best matches chosen by obtaining a score characterising the degree of overlap.



  vector<TrackSegmentCam*> TempTrackSeg[2];
  int Module; bool Cont;

  int id, idm, idp; 
  double TopScore, Score, NClusters, TopClusters;
  int begplane1(0), begplane2(0), endplane1(0), endplane2(0), dplane1, dplane2;
  int NumUTrks(0), NumVTrks(0);
  
  // Store the 2D track segments which roughly overlap in terms of planes
  // and time.
 
  // Loop over all u segments
  for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {
    TrackSegmentCam* Segu = ViewSegBank[0][i];
 
    // If segment is a 2D track   
    if(Segu->GetUID()>2) {

      // For each of these, loop over 2D tracks in V view
      for(unsigned int j=0; j<ViewSegBank[1].size(); ++j) {
        TrackSegmentCam* Segv = ViewSegBank[1][j];
        
        if(Segv->GetUID()>2) {
        
          // Check overlap between u and v segments
          if( (abs(Segv->GetBegPlane()-Segu->GetBegPlane())<10 || abs(Segv->GetEndPlane()-Segu->GetEndPlane())<10 )
              && ( (Segv->GetEndPlane()-Segu->GetBegPlane())>2 && (Segu->GetEndPlane()-Segv->GetBegPlane())>2 )
              && ( (Segv->GetEndTime()-Segu->GetBegTime())>-100 && (Segu->GetEndTime()-Segv->GetBegTime())>-100 ) ) {
            
            TempTrackSeg[0].push_back(Segu);
            TempTrackSeg[1].push_back(Segv);
          }
        }
      }
    }
  }

  // Track combinations we may have missed above:
  // Simple combination
  TrackSegmentCam* USeg=0; TrackSegmentCam* VSeg=0;

  for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {if(ViewSegBank[0][i]->GetUID()>2) {NumUTrks++; USeg=ViewSegBank[0][i];}}
  for(unsigned int i=0; i<ViewSegBank[1].size(); ++i) {if(ViewSegBank[1][i]->GetUID()>2) {NumVTrks++; VSeg=ViewSegBank[1][i];}}

  if(TempTrackSeg[0].size()==0 && TempTrackSeg[1].size()==0 
     && NumUTrks==1 && NumVTrks==1 && USeg && VSeg) {
    TempTrackSeg[0].push_back(USeg); 
    TempTrackSeg[1].push_back(VSeg);
  }


  // Essentially same as above, but if we have also found a junky second track in one view
  if((ModuleType==1 || ModuleType==3) && TempTrackSeg[0].size()==0 && TempTrackSeg[1].size()==0) {
    // First combination
    if(NumUTrks==1 && NumVTrks==2) {
      TrackSegmentCam* VSeg2=0;

      for(unsigned int i=0; i<ViewSegBank[1].size(); ++i) {
        if(ViewSegBank[1][i]->GetUID()>2) {VSeg2=ViewSegBank[1][i]; break;}
      }
      
      if(USeg && VSeg && VSeg2) {
        if(abs(int(USeg->GetEntries())-int(VSeg->GetEntries()))<30 && abs(int(USeg->GetEntries())-int(VSeg2->GetEntries()))>30) {
          TempTrackSeg[0].push_back(USeg); 
          TempTrackSeg[1].push_back(VSeg);
        }
        else if(abs(int(USeg->GetEntries())-int(VSeg->GetEntries()))>30 && abs(int(USeg->GetEntries())-int(VSeg2->GetEntries()))<30) {
          TempTrackSeg[0].push_back(USeg); 
          TempTrackSeg[1].push_back(VSeg2);       
        }
      }      
    }
    // Equivalent combination
    else if(NumUTrks==2 && NumVTrks==1) {
      TrackSegmentCam* USeg2=0;

      for(unsigned int i=0; i<ViewSegBank[0].size(); ++i) {
        if(ViewSegBank[0][i]->GetUID()>2) {USeg2=ViewSegBank[0][i]; break;}
      }
      
      if(VSeg && USeg && USeg2) {
        if(abs(int(VSeg->GetEntries())-int(USeg->GetEntries()))<30 && abs(int(VSeg->GetEntries())-int(USeg2->GetEntries()))>30) {
          TempTrackSeg[0].push_back(USeg); 
          TempTrackSeg[1].push_back(VSeg);
        }
        else if(abs(int(VSeg->GetEntries())-int(USeg->GetEntries()))>30 && abs(int(VSeg->GetEntries())-int(USeg2->GetEntries()))<30) {
          TempTrackSeg[0].push_back(USeg2); 
          TempTrackSeg[1].push_back(VSeg);        
        }
      }      
    }
  }






  // Print out this initial overlap information
  MSG("AlgTrackCamList", Msg::kVerbose) << " ** Guesses at tracks (only showing for UID>2)"<<endl;
  MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
  for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
    TrackSegmentCam* Seg = TempTrackSeg[0][i];
    TrackSegmentCam* SegPartner = TempTrackSeg[1][i];
    if(Seg->GetUID()>2 && SegPartner->GetUID()>2)
      {
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg    " 
                                              << "  begpln: "  << Seg->GetBegPlane()
                                              << "  begtpos: " << Seg->GetBegTPos()
                                              << "  endpln: "  << Seg->GetEndPlane()
                                              << "  endtpos: " << Seg->GetEndTPos() << endl;
        MSG("AlgTrackCamList", Msg::kVerbose) << "  Partner" 
                                              << "  begpln: "  << SegPartner->GetBegPlane()
                                              << "  begtpos: " << SegPartner->GetBegTPos()
                                              << "  endpln: "  << SegPartner->GetEndPlane()
                                              << "  endtpos: " << SegPartner->GetEndTPos() << endl;
        MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
      }
  }




  // Special FD specific section to consider matching 2D tracks across the SM gap
  if(NumModules>1) {

    // Initialise container: Create one vector of empty segment vectors for each SM
    vector< vector<TrackSegmentCam*> > TempMod[2]; 
    for(unsigned int i=0; i<2; ++i) {
      for(unsigned int j=0; j<2; ++j) {vector<TrackSegmentCam*> TempVector; TempMod[i].push_back(TempVector);}
    }

    // Loop over the roughly overlapping 2D tracks we stored above
    for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
      TrackSegmentCam* Segu = TempTrackSeg[0][i];
      TrackSegmentCam* Segv = TempTrackSeg[1][i];
      
      // Check which module the U 2D track ends in
      Module = int((Segu->GetEndPlane())/(PlanesInModule+1));
      
      // Store U and V 2D tracks if they are suitably near the SM gap
      if( Segu->GetBegPlane()<Module*(PlanesInModule+1)+20
          && Segv->GetBegPlane()<Module*(PlanesInModule+1)+20 && Module>0) {
   
        TempMod[1][0].push_back(Segu);
        TempMod[1][1].push_back(Segv);
      }

      if( Segu->GetEndPlane()>(Module+1)*(PlanesInModule+1)-20
          && Segv->GetEndPlane()>(Module+1)*(PlanesInModule+1)-20 && Module<1) {
        
        TempMod[0][0].push_back(Segu);
        TempMod[0][1].push_back(Segv);
      }
    }



    // If there are U and V 2D tracks on either side of the SM gap, look for associations.
    // Also calculate a score based on the degree of overlap we would have if a join was made.
    if(TempMod[0][0].size()>0 && TempMod[1][0].size()>0) {
      Cont=true;

      while(Cont==true) {
        Cont=false;
        idm=-1; idp=-1; TopScore=0.; NClusters=0; TopClusters=0;

        // Loop over the 2D tracks in SM 1
        for(unsigned int i=0; i<TempMod[0][0].size(); ++i) {  
          TrackSegmentCam* Segum = TempMod[0][0][i];
          TrackSegmentCam* Segvm = TempMod[0][1][i];
        
          // For these, loop over the 2D tracks in SM 2
          if( (Segum->GetPartner()==Segvm && Segvm->GetPartner()==Segum)
              || (Segum->GetPartner()==0 && Segvm->GetPartner()==0
                  && Segum->GetUID()>2 && Segvm->GetUID()>2) ) {
            
            for(unsigned int j=0; j<TempMod[1][0].size(); ++j) {
              TrackSegmentCam* Segup = TempMod[1][0][j];
              TrackSegmentCam* Segvp = TempMod[1][1][j];

              // Check the association of the 2D tracks across the gap
              if(Segup->GetUID()>2 && Segvp->GetUID()>2) {
                if(Segum->IsAssoc(Segup) && Segvm->IsAssoc(Segvp)) {
                  
                  if(Segum->GetBegPlane()<=Segvm->GetBegPlane()) {
                    begplane2=Segum->GetBegPlane();
                    begplane1=Segvm->GetBegPlane();
                  }

                  if(Segvm->GetBegPlane()<=Segum->GetBegPlane()) {
                    begplane2=Segvm->GetBegPlane();
                    begplane1=Segum->GetBegPlane();
                  }

                  if(Segup->GetEndPlane()>=Segvp->GetEndPlane()) {
                    endplane2=Segup->GetEndPlane();
                    endplane1=Segvp->GetEndPlane();
                  }

                  if(Segvp->GetEndPlane()>=Segup->GetEndPlane()) {
                    endplane2=Segvp->GetEndPlane();
                    endplane1=Segup->GetEndPlane();
                  }               

                  dplane1=endplane2-begplane2;
                  dplane2=endplane1-begplane1;
                  
                  Score=double(dplane2*dplane2)/double(dplane1);

                  NClusters=int(Segum->GetEntries()+Segvm->GetEntries()+Segup->GetEntries()+Segvp->GetEntries());

                  if(Score>TopScore || (Score==TopScore && NClusters>TopClusters)) {idm=i; idp=j; TopScore=Score; TopClusters=NClusters;}
                }
              }
            }
          }
        }

        // Add 2D tracks together across the SM gap, as a result of the match.
        // Can also reliably set U/V partners. Set UIDs to -1 to avoid considering
        // the old 2D track that is now a subsection of a larger 2D track.
        if((1+idm)>0 && (1+idp)>0) {
          TrackSegmentCam* Segum = TempMod[0][0][idm];
          TrackSegmentCam* Segvm = TempMod[0][1][idm];
          TrackSegmentCam* Segup = TempMod[1][0][idp];
          TrackSegmentCam* Segvp = TempMod[1][1][idp];

          Segup->AddSegment(Segum); Segup->SetPartner(Segvp); 
          Segum->SetUID(-1); Segum->SetPartner(0);
          
          Segvp->AddSegment(Segvm); Segvp->SetPartner(Segup); 
          Segvm->SetUID(-1); Segvm->SetPartner(0);

          Cont=true;
        }
        
      } // End while Cont==true
    }

    // Clean up
    for(unsigned int i=0; i<2; ++i) {TempMod[i].clear();}
  }
  



  // Find other 2D track matches between views.
  Cont=true;
  
  while(Cont==true) {
    Cont=false;
    id=-1; TopScore=0.; NClusters=0; TopClusters=0;

    // Loop over the roughly overlapping 2D tracks we stored above
    for(unsigned int i=0; i<TempTrackSeg[0].size(); ++i) {
      TrackSegmentCam* Segu = TempTrackSeg[0][i];
      TrackSegmentCam* Segv = TempTrackSeg[1][i];

      // If have already established U/V partners, but not made the 3D track
      // do this now. TempTrack separately stores the U and V components of the
      // 3D track, which are currently still track segments, identified by UID==4.
      if(Segu->GetPartner()==Segv && Segv->GetPartner()==Segu
         && Segu->GetUID()<4 && Segv->GetUID()<4) {

        TempTrack[0].push_back(Segu);
        TempTrack[1].push_back(Segv);

        Segu->SetUID(4);
        Segv->SetUID(4);
      }

      // Otherwise, find the 2D tracks which have the best overlap
      else {
        if(Segu->GetPartner()==0 && Segv->GetPartner()==0
           && Segu->GetUID()>2 && Segv->GetUID()>2) {
          
          if(Segu->GetBegPlane()<=Segv->GetBegPlane()) {
            begplane2=Segu->GetBegPlane();
            begplane1=Segv->GetBegPlane();
          }

          if(Segv->GetBegPlane()<=Segu->GetBegPlane()) {
            begplane2=Segv->GetBegPlane();
            begplane1=Segu->GetBegPlane();
          }

          if(Segu->GetEndPlane()>=Segv->GetEndPlane()) {
            endplane2=Segu->GetEndPlane();
            endplane1=Segv->GetEndPlane();
          }

          if(Segv->GetEndPlane()>=Segu->GetEndPlane()) {
            endplane2=Segv->GetEndPlane();
            endplane1=Segu->GetEndPlane();
          }

          dplane1=endplane2-begplane2;
          dplane2=endplane1-begplane1;
                  
          Score=double(dplane2*dplane2)/double(dplane1);

          NClusters=int(Segu->GetEntries()+Segv->GetEntries());

          MSG("AlgTrackCamList", Msg::kVerbose) << "NClusters " << NClusters << " TopClusters "<< TopClusters 
                                                << " Score " << Score << " TopScore " << TopScore << " id " << i << endl;
          
          //      if(Score>TopScore || (Score==TopScore && NClusters>TopClusters)) {id=i; TopScore=Score; TopClusters=NClusters;}
          if((Score>TopScore && NClusters>TopClusters-4)
             || (Score>0.4*TopScore && NClusters>TopClusters+4) 
             || (Score==TopScore && NClusters>TopClusters)) {id=i; TopScore=Score; TopClusters=NClusters;}
        }
      }
    }

    // Set the U/V partners for the 2D tracks having the best overlap
    if((1+id)>0) {
      MSG("AlgTrackCamList", Msg::kVerbose) << "Final id " << id << endl;

      TrackSegmentCam* Segu = TempTrackSeg[0][id];
      TrackSegmentCam* Segv = TempTrackSeg[1][id];
      
      Segu->SetPartner(Segv); 
      Segv->SetPartner(Segu);
      Cont=true;
    }

  } // End while Cont==true



  // Set the cluster flags for those clusters in a 3D track, and print out results
  MSG("AlgTrackCamList", Msg::kVerbose) << " ** Tracks Chosen ** " << endl;
  
  for(int View=0; View<2; ++View) {
    MSG("AlgTrackCamList", Msg::kVerbose) << "View: "<<View<<endl;
    MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;
    
    for(unsigned int i=0; i<TempTrack[View].size(); ++i) {
      TrackSegmentCam* Seg = TempTrack[View][i];
      TrackSegmentCam* SegPartner = Seg->GetPartner();

      MSG("AlgTrackCamList", Msg::kVerbose) << "  Seg    " 
                                            << "  begpln: "  << Seg->GetBegPlane()
                                            << "  begtpos: " << Seg->GetBegTPos()
                                            << "  endpln: "  << Seg->GetEndPlane()
                                            << "  endtpos: " << Seg->GetEndTPos() << endl;
      MSG("AlgTrackCamList", Msg::kVerbose) << "  Partner" 
                                            << "  begpln: "  << SegPartner->GetBegPlane()
                                            << "  begtpos: " << SegPartner->GetBegTPos()
                                            << "  endpln: "  << SegPartner->GetEndPlane()
                                            << "  endtpos: " << SegPartner->GetEndTPos() << endl;
      MSG("AlgTrackCamList", Msg::kVerbose) << "================================" << endl;

      for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
        Seg->GetCluster(j)->SetTrkFlag(4);      
      }
    }
  }

  return;
}




void AlgTrackCamList::Form3DTracks(CandSliceHandle* slice)
{
  // The 3D track segments are now stored in the TempTracks container, with
  // U and V segments stored separately. In this method we look at the consituent
  // clusters and hits in these segments, performing fits to choose the final
  // strips for the track.


  MSG("AlgTrackCamList", Msg::kVerbose) << " *** Formation of 3D tracks *** " << endl;

  int ktemp, View;
  double swx, swy, swx2, swxy, sw, ProjectedTPos, TPosWindow(0.);
  double w, x, y, m, c;
  int UStrips, VStrips, TrackStrips, TrackPlanes;

  int begPlaneView[2], endPlaneView[2];
  int begplane1, begplane2, endplane1, endplane2, maxplane, Plane;
  double dScore, Score; int id;

  int ShwOrTrkList[500];
  int ViewList[500];
  vector<HitCam*> HitList[500];
  vector<TrackCam*> FinalTrackTempHolder[2];
  for(int k=0; k<500; ++k) {ShwOrTrkList[k]=-1; ViewList[k]=-1;}



  // Main loop is over the pairs of track segments, stored above
  for(unsigned int i=0; i<TempTrack[0].size(); ++i) {
    MSG("AlgTrackCamList", Msg::kVerbose) << " making track " << i << endl;
    
    
    // Find the beginning and end planes for each of this pair of segments.
    // Order them so that we have:
    // begplane2, begplane 1, endplane 1, endplane 2.
    begplane2=-1; endplane2=-1; begplane1=-1; endplane1=-1;
    
    for(View=0; View<2; ++View) {
      TrackSegmentCam* seg = TempTrack[View][i];
      if(begplane2<0 || seg->GetBegPlane()<begplane2) { begplane1=begplane2; begplane2=seg->GetBegPlane(); }
      else { begplane1=seg->GetBegPlane(); }
      
      if(endplane2<0 || seg->GetEndPlane()>endplane2) { endplane1=endplane2; endplane2=seg->GetEndPlane(); }
      else { endplane1=seg->GetEndPlane(); }

      begPlaneView[View]=-1; endPlaneView[View]=-1;
    }

    maxplane=1+endplane2;
    
    MSG("AlgTrackCamList", Msg::kVerbose) << begplane2 << " (" << begplane1 << ") -> (" << endplane1 << ") " << endplane2 << endl;
 

    // Unpack the clusters from the segments.
    // For each track, should only have one cluster on a plane.
    for(View=0; View<2; ++View) {
      TrackSegmentCam* Seg = TempTrack[View][i];
      
      for(unsigned int j=0; j<Seg->GetEntries(); ++j) {
        ClusterCam* Clust = Seg->GetCluster(j);

        Plane=Clust->GetPlane();
        if(Plane>=begplane2 && Plane<maxplane) {ClusterList[Plane].push_back(Clust);}
      }
    }
    
    
    
    // Trim the segments. Check we haven't tracked much further
    // in one view than the other, and adjust as necessary.

    // Maximum discrepancy allowed at each end is 3 planes. As only
    // clusters on planes greater than begplane1 are considered.
    int PrevPln, NextPln;
    bool LeavesAtEdge=false;

    // Check hasn't just left edge of detector in one view.
    if(ModuleType==1 && (ClusterList[begplane1][0]->GetEndTPos()>3.9 
                         || ClusterList[begplane1][0]->GetBegTPos()<-3.9)) {begplane1=begplane2+1; LeavesAtEdge=true;}
    if(ModuleType==1 && (ClusterList[endplane1][0]->GetEndTPos()>3.9 
                         || ClusterList[endplane1][0]->GetBegTPos()<-3.9)) {endplane1=endplane2-1; LeavesAtEdge=true;}

    if(!LeavesAtEdge) {
      // For the beginning
      // Jump back 4 planes before begplane1 (unless this takes us past begplane2)
      begplane1=begplane1-4; PrevPln=-1;
      Plane=begplane1; if(Plane<begplane2) {Plane=begplane2;}
    
      // Working forwards from this plane, find the first plane that has
      // a cluster.
      while(Plane<maxplane) {
        if(ClusterList[Plane].size()>0) {PrevPln=Plane; break;}
        else {Plane++;}
      }
    
      // Look at the next 2 planes in the same view as this cluster
      // and see if they have clusters too. Try and veto geometries such as
      //
      //                bp1
      // u       x   x   o   o
      // v         o   x   *   o
      //          bp2      * == shw-like, or gap
      //
      // In this case, begplane1 is set to begplane 2, so the lone hit at 
      // begplane 2 is later vetoed.
      if(PrevPln>-1 && PrevPln+4<maxplane && ClusterList[PrevPln+2].size()==0) {
        if(ClusterList[PrevPln+4].size()>0) {
          ClusterCam* Clust = ClusterList[PrevPln+4][0];
        
          if(Clust->GetTrkPlnFlag()>0 && Clust->GetShwPlnFlag()==0) {
            PrevPln=begplane1;
          }
        }
        begplane1=PrevPln;
      }
    
    
      // Similarly, for the end 
      // Jump forward 4 planes after endplane1 (unless this takes us past endplane2)
      endplane1=endplane1+4; NextPln=-1;
      Plane=endplane1; if(Plane>maxplane-1) {Plane=maxplane-1;}
    
      // Working backwards from this plane, find the first plane that has
      // a cluster.
      while(Plane>begplane2-1) {
        if(ClusterList[Plane].size()>0) {NextPln=Plane; break;}
        else {Plane--;}
      }
    
      // Look at the next 2 planes in the same view as this cluster
      // and see if they have clusters too, as above.
      if(NextPln>-1 && NextPln-4>(begplane2-1) && ClusterList[NextPln-2].size()==0) {
        if(ClusterList[NextPln-4].size()>0) {
          ClusterCam* Clust = ClusterList[NextPln-4][0];
        
          if(Clust->GetTrkPlnFlag()>0 && Clust->GetShwPlnFlag()==0) {
            NextPln=endplane1;
          }
        }
        endplane1=NextPln;
      }
    }

    

    // Count the strips and record planeviews, number of occupied planes, etc.
    UStrips=0; VStrips=0; TrackStrips=0;

    // Look at the region over which the track is now defined
    for(int k=begplane2; k<maxplane; ++k) {  

      if( //  (k>=begplane1 && k<=endplane1) &&
          ClusterList[k].size()>0 ) {

        ClusterCam* Clust = ClusterList[k][0];
        
        View=Clust->GetPlaneView();
        ViewList[k]=View;
        ShwOrTrkList[k]=1;
        
        // Record first and last planes in each view.
        if(begPlaneView[View]<0 || k<begPlaneView[View]) {begPlaneView[View]=k;}
        if(endPlaneView[View]<0 || k>endPlaneView[View]) {endPlaneView[View]=k;}
        
        if(View==0) {UStrips++;}
        if(View==1) {VStrips++;}
        
        if(Clust->GetTrkPlnFlag()>0) {TrackStrips++;}
      }
    }
    
    
    MSG("AlgTrackCamList", Msg::kVerbose) 
      << " U: (" << begPlaneView[0] << "," << endPlaneView[0] << ") " 
      << " V: (" << begPlaneView[1] << "," << endPlaneView[1] << ") " 
      << " UStrips " << UStrips << " VStrips " << VStrips << " TrackStrips " << TrackStrips << endl;
    
    

    
    // If we have at least 3 strips in each view and at least 2 clusters
    // on track-like planes, proceed with the fits.
    if(UStrips>2 && VStrips>2 && (CambridgeAnalysis==false || TrackStrips>2) ) {


      // FLAG PLANES which we think are track-like
      for(View=0; View<2; ++View) {
        for(Plane=begplane2; Plane<maxplane; ++Plane) {

          if(ViewList[Plane]==View && ShwOrTrkList[Plane]>0) {
            ClusterCam* Clust0 = ClusterList[Plane][0];

            // Always set flag to 3 for the +/- 5 plane ND track segments
            if(Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane]=3;}


            // Find next planes, in the same view, where we have a cluster, both forwards
            // and backwards.
            ktemp=0; PrevPln=-1;
            while(Plane-ktemp>begplane2) {
              ktemp++;
              if(ViewList[Plane-ktemp]==View) {PrevPln=ktemp; break;}
            }

            ktemp=0; NextPln=-1;
            while(Plane+ktemp<maxplane-1) {
              ktemp++;
              if(ViewList[Plane+ktemp]==View) {NextPln=ktemp; break;}
            }


            // If there are planes forwards and backwards in the same view...
            if(PrevPln>0 && NextPln>0) {
              ClusterCam* NextClust = ClusterList[Plane+NextPln][0];
              ClusterCam* PrevClust = ClusterList[Plane-PrevPln][0];
              
              // If either both clusters are within range for this view, or they are the
              // next plane along and the current plane is a track plane.
              if( (PrevClust->GetPlane()>begPlaneView[View]
                   || (Clust0->GetPlane()-PrevClust->GetPlane()<3 && Clust0->GetTrkPlnFlag()>0))
                  && (NextClust->GetPlane()<endPlaneView[View]
                      || (NextClust->GetPlane()-Clust0->GetPlane()<3 && Clust0->GetTrkPlnFlag()>0) )) {


                // If the clusters have track-like association, identify plane as
                // being track like, by setting entry in ShwOrTrkList to 3.
                if(Clust0->IsTrkAssoc(PrevClust,NextClust)>1) {
                  ShwOrTrkList[Plane]=3;

                  // If there is a good overlap of strip numbers between the three clusters
                  if( (NextClust->GetBegTPos()-Clust0->GetEndTPos()>-1.1*StripWidth 
                       && Clust0->GetBegTPos()-PrevClust->GetEndTPos()>-1.1*StripWidth)
                      || (NextClust->GetEndTPos()-Clust0->GetBegTPos()<1.1*StripWidth
                          && Clust0->GetEndTPos()-PrevClust->GetBegTPos()<1.1*StripWidth) ) {

                    // If forwards and backwards planes are adjacent in view,
                    // identify these planes as track-like too
                    if( Clust0->GetPlane()-PrevClust->GetPlane()<3 || Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane-PrevPln]=3;}
                    if( NextClust->GetPlane()-Clust0->GetPlane()<3 || Clust0->GetNDFlag()==2) {ShwOrTrkList[Plane+NextPln]=3;}
                  }
                }

                // Otherwise, if cluster is small, can identify its plane as track-like.
                // Don't set anything for forwards and backwards planes.
                else {
                  if(Clust0->GetEndTPos()-Clust0->GetBegTPos()<2.1*StripWidth 
                     && Clust0->GetTrkPlnFlag()>0) {ShwOrTrkList[Plane]=3;}
                }
              }
            }

          }
        } // End loop over planes
      } // End loop over views


      MSG("AlgTrackCamList", Msg::kVerbose) << " *** SORT TRACK HITS *** " << endl;
      
      for(Plane=begplane2;Plane<maxplane; ++Plane) {
        if(ViewList[Plane]>-1) {
          ClusterCam* Clust = ClusterList[Plane][0];
          MSG("AlgTrackCamList", Msg::kVerbose) << Clust->GetPlane() 
                                                << " " << Clust->GetBegTPos() << " " << Clust->GetEndTPos()
                                                << " " << ShwOrTrkList[Plane] << endl;
        }
      }



      // Perform GLOBAL LINEAR FIT
      int ShowerStart, ShowerEnd, FitStart, FitEnd, PrevTrkPln, NextTrkPln;

      for(View=0; View<2; ++View) {
        for(Plane=begplane2; Plane<maxplane; ++Plane) {


          // Currently ShwOrTrkList flag==1 means the plane was part of the 3D track.
          // ShwOrTrkList flag==3 means that we also flagged the plane as being track-like, above.
          //
          // Here, we look at the planes we haven't flagged as track-like, find track-like planes
          // on either side and carry out a linear fit through this 'shower-like' region.
          // 
          // Using the fit, clusters in this shower-like region are identified, then hits.
          if(ViewList[Plane]==View && ShwOrTrkList[Plane]<2) {
            MSG("AlgTrackCamList", Msg::kVerbose) << " *** GLOBAL FIT ***   Plane=" << Plane << endl;

            // Look for next planes in the view that we have flagged as track-like, 
            // forwards and backwards.
            ktemp=0; NextTrkPln=-1;
            while(Plane+ktemp<endplane1 && Plane+ktemp<maxplane-1) {
              ktemp++;
              if(ViewList[Plane+ktemp]==View && ShwOrTrkList[Plane+ktemp]>2) {NextTrkPln=ktemp; break;}
            }
            
            ktemp=0; PrevTrkPln=-1;
            while(Plane-ktemp>begplane1 && Plane-ktemp>begplane2) {
              ktemp++;
              if(ViewList[Plane-ktemp]==View && ShwOrTrkList[Plane-ktemp]>2) {PrevTrkPln=ktemp; break;}
            }

            
            // Find the track-like planes which are just outside of the shower-like region ('ShowerStart' 
            // and 'ShowerEnd' planes away), as well as the planes which are up to 7 planes either side 
            // of the shower-like region ('FitStart' and 'FitEnd' planes away). These are the 
            // track used in the linear fit through the shower.
            if(PrevTrkPln>0) {FitStart=-PrevTrkPln-7; ShowerStart=-PrevTrkPln;} 
            else {FitStart=begplane1-Plane+1; ShowerStart=FitStart;}

            if(Plane+FitStart<begplane2) {FitStart=-Plane+begplane2;}
            if(Plane+ShowerStart<begplane2) {ShowerStart=-Plane+begplane2;}

            if(NextTrkPln>0) {FitEnd=1+NextTrkPln+7; ShowerEnd=1+NextTrkPln;}
            else {FitEnd=endplane1-Plane; ShowerEnd=FitEnd;}

            if(Plane+FitEnd>maxplane) {FitEnd=maxplane-Plane;}
            if(Plane+ShowerEnd>maxplane) {ShowerEnd=maxplane-Plane;}


            ClusterCam* ClustSeed = ClusterList[Plane][0];

            m=0.; c=0.5*(ClustSeed->GetBegTPos()+ClustSeed->GetEndTPos());
            swx=0.; swy=0.; swx2=0.; swxy=0.; sw=0.;


            // Carry out the linear fit through the shower, obtaining gradient and intercept.
            for(int k=FitStart; k<FitEnd; ++k) {
              if(ViewList[Plane+k]==View) {
                ClusterCam* ClustTemp = ClusterList[Plane+k][0];
                
                for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
                  HitCam* HitTemp = ClustTemp->GetHit(l);

                  x=HitTemp->GetZPos();
                  y=HitTemp->GetTPos();
                  w=HitTemp->GetCharge()/ClustTemp->GetCharge();
                  
                  if(ShwOrTrkList[Plane+k]>2) {w=5.*w;}
                  else {w=0.5*w;}
                  
                  swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
                }
              }
            }

            if(sw>0. && (sw*swx2-swx*swx)!=0.) {
              m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
              c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
            }


            // Now loop over the planes in the shower-like region
            for(int k=ShowerStart; k<ShowerEnd; ++k) {
              if(ViewList[Plane+k]==View && ShwOrTrkList[Plane+k]<2) {
                // ***First cluster in ClusterList is that from segment in 3D track***
                ClusterCam* ClustTemp = ClusterList[Plane+k][0];
                
                // Projected strip number in shower
                ProjectedTPos = m*ClustTemp->GetZPos()+c;

                // Calculate window around projected strip number
                if(m>=0.) {TPosWindow = (0.6*StripWidth)+0.02*m;}
                if(m<=0.) {TPosWindow = (0.6*StripWidth)-0.02*m;}

                MSG("AlgTrackCamList", Msg::kVerbose) << " Plane=" << Plane+k << " ProjectedTPos=" << ProjectedTPos 
                                                      << " TPosWindow=" << TPosWindow << " (" << ClustTemp->GetBegTPos() 
                                                      << " " << ClustTemp->GetEndTPos() << ") " << endl;


                // If  1. we have found track-like planes on either side of the shower-like region,
                // or  2. we haven't found any of these track-like planes (entire track in this view must be shower),
                // or  3. current cluster overlaps quite nicely with projected strip number from the fit,
                //
                // then proceed with trying to find the best hits from within the cluster.
                if( (PrevTrkPln>0 && NextTrkPln>0) || (PrevTrkPln<0 && NextTrkPln<0)
                    || (ClustTemp->GetEndTPos()>ProjectedTPos-(1.6*StripWidth) && ClustTemp->GetBegTPos()<ProjectedTPos+(1.6*StripWidth)) ) {

                  id=-1; dScore=0.; Score=999.9;

                  // Loop over hits in cluster and get a 'score' for each. Score is based
                  // on distance from projected tpos, and we want to minimise this.
                  for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
                    HitCam* HitTemp = ClustTemp->GetHit(l);
                    
                    dScore=HitTemp->GetTPos()-ProjectedTPos;
                    if(dScore<0) {dScore=-dScore;}

                    if(dScore<Score) {id=l; Score=dScore;}
                  }

        
                  // If we have found a 'best hit'
                  if((1+id)>0) {
                    HitCam* Hit = ClustTemp->GetHit(id);
                    
                    // Create new cluster and store it. 
                    // ***Last cluster in ClusterList will be that found by the fit.***
                    ClusterCam* ClustNew = new ClusterCam(Hit);
                    ClusterList[Plane+k].push_back(ClustNew);
                    
                    ShwOrTrkList[Plane+k]=2;
                    
                    // If any other hits in the cluster found by the fit are also good, 
                    // add these to the new cluster too.
                    for(unsigned int l=0; l<ClustTemp->GetEntries(); ++l) {
                      HitCam* HitTemp = ClustTemp->GetHit(l);
                        
                      if( (HitTemp->GetTPos()-Hit->GetTPos()>0
                           && HitTemp->GetTPos()-Hit->GetTPos()<TPosWindow)
                          || (HitTemp->GetTPos()-Hit->GetTPos()<0
                              && HitTemp->GetTPos()-Hit->GetTPos()>-TPosWindow) ) {
                          
                        ClustNew->AddHit(HitTemp);
                      }
                    }

                    for(unsigned int l=0; l<ClustNew->GetEntries(); ++l) {
                      MSG("AlgTrackCamList", Msg::kVerbose) << "   ...global add hit=" << ClustNew->GetHit(l)->GetTPos() << " (" << Plane+k << ") " << endl;
                    }
                  }

                }
              }
            }
          }
        }  // End loop over planes
      } // End loop over views



      // Now, Perform LOCAL LINEAR FITS
      for(View=0; View<2; ++View) {
        for(Plane=begplane2; Plane<maxplane; ++Plane) {

          // If plane has been declared track-like
          if(ViewList[Plane]==View && ShwOrTrkList[Plane]>1) {
            MSG("AlgTrackCamList", Msg::kVerbose) << " *** FIT (2) ***   Plane=" << Plane << endl;

            // Get last cluster from ClusterList. This will have been flagged immediately as track-like,
            // or have been stored as a result of a linear fit through a shower-like region.
            // Using the cluster, calculate a projected strip number and window for use in the local fit.
            ClusterCam* clr = ClusterList[Plane].back();
            ProjectedTPos=0.5*(clr->GetBegTPos()+clr->GetEndTPos()); 
            TPosWindow=0.6*StripWidth;

            // Loop over the hits in this cluster
            const unsigned int nhits = clr->GetEntries();
            if(nhits>1){
              for(unsigned int k=0; k<nhits; ++k) {clr->GetHit(k)->SetTrkFlag(1);}
              
              swx=0.0; swy=0.0; swx2=0.0; swxy=0.0; sw=0.0;

              // For these, loop over the hits in nearby clusters (within 4 planes of the current plane)
              // and carry out the linear fit
              for(int k=-4; k<5; ++k){
                if( (Plane+k)>=begplane2 && (Plane+k)<(maxplane) && ViewList[Plane+k]==View && ShwOrTrkList[Plane+k]>1 ) {
                  ClusterCam* clrtmp = ClusterList[Plane+k].back(); // Again get last entry 

                  for(unsigned int l=0; l<clrtmp->GetEntries(); ++l){
                    HitCam* hittmp = clrtmp->GetHit(l);
                    x=hittmp->GetZPos(); 
                    y=hittmp->GetTPos(); 
                    w=hittmp->GetCharge()/clrtmp->GetCharge();

                    swx+=w*x; swy+=w*y; swx2+=w*x*x; swxy+=w*x*y; sw+=w;
                  }
                }
              }

              if(sw>1.0 && (sw*swx2-swx*swx)!=0.){
                m=(sw*swxy-swx*swy)/(sw*swx2-swx*swx);
                c=(swy*swx2-swx*swxy)/(sw*swx2-swx*swx);
                ProjectedTPos=m*clr->GetZPos()+c; 
                if(m>=0.0) TPosWindow = (0.6*StripWidth)+0.02*m; 
                if(m<=0) TPosWindow = (0.6*StripWidth)-0.02*m; 
              }
            }
            
            MSG("AlgTrackCamList", Msg::kVerbose) << "  local plane=" << Plane << " ProjectedTPos=" << ProjectedTPos 
                                                  << " TPosWindow=" << TPosWindow << endl;


            // Using the results of the fit, obtain a score for each hit. Again, this is
            // based on distance from the projected strip number. We want to pick the hit which
            // minimises this distance.
            id=-1; dScore=0.; Score=999.;
            
            for(unsigned int k=0; k<nhits; ++k){
              HitCam* hit = clr->GetHit(k);  
              dScore=hit->GetTPos()-ProjectedTPos; if(dScore<0) dScore=-dScore;
              
              if(dScore<Score){id=k; Score=dScore;}
            }
            
            // If we have found a good hit, store it and flag the plane as part of the 
            // final track
            if(1+id>0){  
              HitCam* hit = clr->GetHit(id); 
              HitList[Plane].push_back(hit); ShwOrTrkList[Plane]=4;

              // Include any other hits which also match projection within the window
              for(unsigned int k=0; k<nhits; ++k){
                HitCam* hittmp = clr->GetHit(k);
                if( ( hittmp->GetTPos()-hit->GetTPos()>0 && hittmp->GetTPos()-hit->GetTPos()<TPosWindow)
                    || ( hittmp->GetTPos()-hit->GetTPos()<0 && hittmp->GetTPos()-hit->GetTPos()>-TPosWindow) ){
                  HitList[Plane].push_back(hittmp);
                }
              }
            }

            for(unsigned int l=0; l<HitList[Plane].size(); ++l){
              HitCam* hittmp = HitList[Plane][l];
              MSG("AlgTrackCamList", Msg::kVerbose) << "   ...local add hit=" << hittmp->GetTPos() << " (" << Plane << ") " << endl;
            }

          } // End demand that plane has been flagged as track-like
        } // End loop over planes
      } // End loop over views
      


      // *** FINALLY, FORM THE TRACKS ***
      // Count the number of planes flagged as being part of final track
      TrackPlanes=0;
      
      for(Plane=begplane2; Plane<maxplane; ++Plane) {
        if(ShwOrTrkList[Plane]>3) {TrackPlanes++;}
      }

      MSG("AlgTrackCamList", Msg::kVerbose) << " TrackPlanes " << TrackPlanes << endl;

      if(CambridgeAnalysis==false || TrackPlanes>5) {
        for(View=0; View<2; ++View) {
          
          TrackCam* Trk = new TrackCam(slice);
          FinalTrackTempHolder[View].push_back(Trk);
          
          for(Plane=begplane2; Plane<maxplane; ++Plane) {
            if(ViewList[Plane]==View && ShwOrTrkList[Plane]>3) {
              
              for(unsigned int k=0; k<HitList[Plane].size(); ++k) {
                HitList[Plane][k]->SetTrkFlag(2);
                Trk->AddHit(HitList[Plane][k]);
              }
            }
          }
        }
        
        TrackCam* Trku = FinalTrackTempHolder[0].back();
        TrackCam* Trkv = FinalTrackTempHolder[1].back();
        
        Trku->SetPartner(Trkv); 
        Trkv->SetPartner(Trku);
        
        for(unsigned int p=0; p<Trku->GetEntries(); ++p) {HitsInTracks[0].push_back(Trku->GetHit(p));}
        for(unsigned int p=0; p<Trkv->GetEntries(); ++p) {HitsInTracks[1].push_back(Trkv->GetHit(p));}
      }
      
    } // End U/V/TrackStrips check

    for(int k=0; k<500; ++k) {
      ShwOrTrkList[k]=-1; 
      ViewList[k]=-1;
      ClusterList[k].clear();
      HitList[k].clear();
    }

  } // End loop over temporary tracks   
 
 

  // Check tracks don't contain same hits as a previous track
  bool GoodTrack; bool ThreePlanes[2];
  int DuplicateStrips[2]; int Planes[2]; int ThisPlane;
  vector<HitCam*> HitsInPrevTracks[2];

  for(unsigned int p=0; p<FinalTrackTempHolder[0].size(); ++p) {

    TrackCam* Trku = FinalTrackTempHolder[0][p];
    TrackCam* Trkv = FinalTrackTempHolder[1][p];

    if(ModuleType==1) {LookForHitsAcrossGap(Trku); LookForHitsAcrossGap(Trkv);}

    ExtendTrack(Trku); ExtendTrack(Trkv);
    FillGapsInTrack(Trku); FillGapsInTrack(Trkv);


    GoodTrack=false; 

    if(Trku->GetEntries()>2 && Trkv->GetEntries()>2) {
      DuplicateStrips[0]=0; DuplicateStrips[1]=0;

      // For u
      Planes[0]=-999; Planes[1]=-999; ThreePlanes[0]=false;

      for(unsigned int p=0; p<Trku->GetEntries(); ++p) {
        HitCam* TrackHit = Trku->GetHit(p);

        for(unsigned int q=0; q<HitsInPrevTracks[0].size(); ++q) {
          if(TrackHit==HitsInPrevTracks[0][q]) {DuplicateStrips[0]++; break;}
        }

        ThisPlane=TrackHit->GetPlane();
        if(Planes[0]<0) {Planes[0]=ThisPlane;}
        else if(Planes[1]<0 && Planes[0]!=ThisPlane) {Planes[1]=ThisPlane;}
        else if(Planes[0]!=ThisPlane && Planes[1]!=ThisPlane) {ThreePlanes[0]=true;}

        if(DuplicateStrips[0]>2) {break;}
      }


      // For v
      Planes[0]=-999; Planes[1]=-999; ThreePlanes[1]=false;

      for(unsigned int p=0; p<Trkv->GetEntries(); ++p) {
        HitCam* TrackHit = Trkv->GetHit(p);

        for(unsigned int q=0; q<HitsInPrevTracks[1].size(); ++q) {
          if(TrackHit==HitsInPrevTracks[1][q]) {DuplicateStrips[1]++; break;}
        }

        ThisPlane=TrackHit->GetPlane();
        if(Planes[0]<0) {Planes[0]=ThisPlane;}
        else if(Planes[1]<0 && Planes[0]!=ThisPlane) {Planes[1]=ThisPlane;}
        else if(Planes[0]!=ThisPlane && Planes[1]!=ThisPlane) {ThreePlanes[1]=true;}

        if(DuplicateStrips[1]>2) {break;}
      }

      if(ThreePlanes[0] && ThreePlanes[1] && DuplicateStrips[0]<3 && DuplicateStrips[1]<3) {GoodTrack=true;}    
    }


    if(GoodTrack) {
      FinalTrackBank[0].push_back(Trku);
      FinalTrackBank[1].push_back(Trkv);

      for(unsigned int q=0; q<Trku->GetEntries(); ++q) {HitsInPrevTracks[0].push_back(Trku->GetHit(q));}
      for(unsigned int q=0; q<Trkv->GetEntries(); ++q) {HitsInPrevTracks[1].push_back(Trkv->GetHit(q));}

      if(ModuleType==2) {MatchUV(Trku,Trkv);}
    }
  }

  FinalTrackTempHolder[0].clear(); FinalTrackTempHolder[1].clear();
  HitsInPrevTracks[0].clear(); HitsInPrevTracks[1].clear();



  // Print out list of final tracks           
  MSG("AlgTrackCamList", Msg::kVerbose) << " *** LIST OF TRACKS *** " << endl;
  for(View=0; View<2; ++View) {
    for(unsigned int i=0; i<FinalTrackBank[View].size(); ++i) {
      TrackCam* TrkTemp = FinalTrackBank[View][i];
      
      MSG("AlgTrackCamList", Msg::kVerbose) 
        << "  " << View << " " << i
        << "  " << TrkTemp->GetBegPlane() << "  " << TrkTemp->GetEndPlane() << endl;
    }
  }

  return;
}



void AlgTrackCamList::ExtendTrack(TrackCam* Trk)
{
  MSG("AlgTrackCamList", Msg::kVerbose) << " Attempting to extend track at beginning and end " << endl;
  
  int PlanesSinceAdded, CurrentPlane(0), Increment(0);
  unsigned int nhits;
  double TrkDir(0.), TrkTPos(0.), TrkZPos(0.), PredictedTPos, extrem1, extrem2;
  double StripCentre, CurrentZPos, MinDistanceToHit, StripCharge;

  double Tolerance; int MaxPlanes;
  int PlaneView=Trk->GetPlaneView();

  bool ConsiderHit;

  if(ModuleType==1 || ModuleType==3) {MaxPlanes=8; Tolerance=2.*StripWidth;}
  else {MaxPlanes=40; Tolerance=4.*StripWidth;}

  // Extend beginning first, then the end.
  for(int direction=0; direction<2; ++direction) {
  
    PlanesSinceAdded=1; 
    if(direction==0) {CurrentPlane=Trk->GetBegPlane()-1; Increment=-1;}
    else if(direction==1) {CurrentPlane=Trk->GetEndPlane()+1; Increment=1;}

    
    while(PlanesSinceAdded<=MaxPlanes && CurrentPlane>=0 && CurrentPlane<=PlanesInModule) { 
      // Don't do anything if there are no hits, or we aren't in the same view.
      nhits=AllHitBank[CurrentPlane].size();
      if(nhits>0 && AllHitBank[CurrentPlane][0]->GetPlaneView()==PlaneView) {

        if(direction==0) {
          TrkDir=Trk->GetBegDir(); 
          TrkTPos=Trk->GetBegTPos(); 
          TrkZPos=Trk->GetBegZPos(); 
        }
        else if(direction==1) {
          TrkDir=Trk->GetEndDir(); 
          TrkTPos=Trk->GetEndTPos(); 
          TrkZPos=Trk->GetEndZPos(); 
        }

        HitCam* CurrentHit = 0;
        CurrentZPos = AllHitBank[CurrentPlane][0]->GetZPos();

        MinDistanceToHit = Tolerance;
        
        // Loop over the hits on the plane
        for(unsigned int i=0; i<nhits; ++i) {

          // Don't add a hit that is already in another track in the slice
          ConsiderHit=true;
          for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
            if(AllHitBank[CurrentPlane][i]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
          }
          if(ConsiderHit==false) {continue;}


          // Find the most likely hit to add
          PredictedTPos = TrkTPos + TrkDir*(CurrentZPos-TrkZPos);
          extrem1 = PredictedTPos + (0.055*TrkDir);
          extrem2 = PredictedTPos - (0.055*TrkDir);
          
          StripCentre=AllHitBank[CurrentPlane][i]->GetTPos();
          StripCharge=AllHitBank[CurrentPlane][i]->GetCharge();
          
          MSG("AlgTrackCamList", Msg::kVerbose) << "Prediction " << direction << " " << PredictedTPos 
                                                << " Actual " << StripCentre   << " CurrentZPos " << CurrentZPos << " trkdir " << TrkDir << endl;
          
          if(fabs(StripCentre-PredictedTPos)<MinDistanceToHit 
             && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-PredictedTPos)<StripWidth && PlanesSinceAdded<20)) ) {
            MinDistanceToHit = fabs(StripCentre-PredictedTPos);
            CurrentHit = AllHitBank[CurrentPlane][i];
          }
          else if(fabs(StripCentre-extrem1)<MinDistanceToHit 
                  && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-extrem1)<StripWidth && PlanesSinceAdded<20)) ) {
            MinDistanceToHit = fabs(StripCentre-extrem1);
            CurrentHit = AllHitBank[CurrentPlane][i];
          }
          else if(fabs(StripCentre-extrem2)<MinDistanceToHit 
                  && (ModuleType==1 || ModuleType==3 || StripCharge>2 || (fabs(StripCentre-extrem2)<StripWidth && PlanesSinceAdded<20)) ) {
            MinDistanceToHit = fabs(StripCentre-extrem2);
            CurrentHit = AllHitBank[CurrentPlane][i];
          }
        }
        
        // Add best hit we have found
        if(CurrentHit) {
          Trk->AddHit(CurrentHit); 
          HitsInTracks[PlaneView].push_back(CurrentHit);

          PlanesSinceAdded=0;

          MSG("AlgTrackCamList", Msg::kVerbose) << "Track extended " << direction << " TPos " << CurrentHit->GetTPos() 
                                                << " ZPos " << CurrentHit->GetZPos() << endl;
        }
      }
      else {PlanesSinceAdded++;}
      
      CurrentPlane+=Increment;
    }
  }


  return;
}



void AlgTrackCamList::FillGapsInTrack(TrackCam* Trk)
{
  MSG("AlgTrackCamList", Msg::kVerbose) << " Attempting to fill gaps in track " << endl;


  // First store track hits in plane order
  vector<HitCam*> TrackHits[490];
  int PlaneView=Trk->GetPlaneView();

  for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
    HitCam* Hit = Trk->GetHit(i);
    TrackHits[Hit->GetPlane()].push_back(Hit);
  }

  
  // Find gaps
  int PrevPlaneWithHit=Trk->GetBegPlane();
  int NextPlaneWithHit=Trk->GetEndPlane();

  double PrevTPos, NextTPos, PrevZPos, NextZPos, CurrentZPos;
  double PredictedGradient, PredictedTPos, MinDistanceToHit, StripCharge;
  double Tolerance; bool ConsiderHit;

  if(ModuleType==1 || ModuleType==3) {Tolerance=1.5*StripWidth;}
  else {Tolerance=3.*StripWidth;}

  for(int i=Trk->GetBegPlane()+1; i<Trk->GetEndPlane(); i++) {
    // Don't do anything if there are no hits, or we aren't in the same view.
    if(AllHitBank[i].size()>0 && AllHitBank[i][0]->GetPlaneView()==PlaneView) {
      
      // Find the hit planes on either side of a gap.
      if(TrackHits[i].size()>0) {PrevPlaneWithHit=i; continue;}

      for(int j=i+2; j<=Trk->GetEndPlane(); j++) {
        if(TrackHits[j].size()>0) {NextPlaneWithHit=j; break;}
      }
 
      // Now try and fill the gap. Limit to maximum gap size?
      if(PrevPlaneWithHit!=NextPlaneWithHit) {
        CurrentZPos=AllHitBank[i][0]->GetZPos();

        PrevZPos=TrackHits[PrevPlaneWithHit][0]->GetZPos();
        NextZPos=TrackHits[NextPlaneWithHit][0]->GetZPos();

        PrevTPos=TrackHits[PrevPlaneWithHit][0]->GetTPos();
        NextTPos=TrackHits[NextPlaneWithHit][0]->GetTPos();

        PredictedGradient=(NextTPos-PrevTPos)/(NextZPos-PrevZPos);
        PredictedTPos=PrevTPos+((CurrentZPos-PrevZPos)*PredictedGradient);

        HitCam* CurrentHit = 0;

        MinDistanceToHit=Tolerance;

        for(unsigned int k=0; k<AllHitBank[i].size(); ++k) {
          HitCam* hit = AllHitBank[i][k];


          // Don't add a hit that is already in another track in the slice
          ConsiderHit=true;
          for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
            if(hit==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
          }
          if(ConsiderHit==false) {continue;}
          

          // Find the most likely hit to add
          StripCharge=hit->GetCharge();

          if(fabs(PredictedTPos-hit->GetTPos())<MinDistanceToHit 
             && (ModuleType==1 || ModuleType==3 || StripCharge>2 || fabs(PredictedTPos-hit->GetTPos())<StripWidth)) {

            MinDistanceToHit=fabs(PredictedTPos-hit->GetTPos());
            CurrentHit=hit;
          }
        }
      
        if(CurrentHit) {
          Trk->AddHit(CurrentHit); 
          TrackHits[i].push_back(CurrentHit);
          HitsInTracks[PlaneView].push_back(CurrentHit);

          MSG("AlgTrackCamList", Msg::kVerbose) << "Gap filled at TPos " << CurrentHit->GetTPos()
                                                << " ZPos " << CurrentHit->GetZPos() << endl;
        }
      }
    }
  }
  
  // Clean up
  for(int i=0; i<490; ++i) {TrackHits[i].clear();}
  
  return;
}




void AlgTrackCamList::LookForHitsAcrossGap(TrackCam* Trk)
{
  // Try and add missing hits that are just across the SM gap
  
  int PlanesSinceAdded, CurrentPlane(0), HitsAdded(0);
  double CurrentZPos, TargetTPos;
  double TrkBegZPos, TrkBegTPos, TrkBegDir(0.);
  double TrkEndZPos, TrkEndTPos, TrkEndDir(0.);
  double Tolerance(0.);

  double MinDistanceToHit, DistanceToCurrentHit;
  double centre, extrem1, extrem2;
  unsigned int nhits; 

  double z, t, w, sw, swx, swx2, swy, swyx;
  bool ConsiderHit;
  
  int PlaneView=Trk->GetPlaneView();

  // If we want to look back into SM 1
  if(Trk->GetBegPlane()>248 && Trk->GetBegPlane()<=256) {
              
    if(PlaneView==0) {CurrentPlane=248;} else if(PlaneView==1) {CurrentPlane=247;}
    PlanesSinceAdded=0;
    TrkBegZPos=Trk->GetBegZPos();
    TrkBegTPos=Trk->GetBegTPos();
 

    // Get beginning direction
    z=0.; t=0.; sw=0.; swx=0.; swx2=0.; swy=0.; swyx=0.;
     
    // Include first few hits in track
    for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
      HitCam* hit = Trk->GetHit(i);

      if(hit->GetPlane()<Trk->GetBegPlane()+8) {
        z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
        sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
      }
    }
    // Use result so far to define a tolerance
    if((swx*swx-sw*swx2)!=0) {TrkBegDir = (swx*swy-sw*swyx)/(swx*swx-sw*swx2); Tolerance = 0.2 + fabs(TrkBegDir);}

              
    // Look at some of the relevant hits across the gap
    for(int i=CurrentPlane; i>CurrentPlane-8; i-=2) {
      for(unsigned int j=0; j<AllHitBank[i].size(); ++j) {
        HitCam* hit = AllHitBank[i][j];
          
        z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();

        // If hit is within tolerance region, include with weight 0.5
        if(fabs(t-(TrkBegTPos+TrkBegDir*(z-TrkBegZPos)))<Tolerance) {
          sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
        }
      }
    }

    // Get results of total fit
    if((swx*swx-sw*swx2)!=0) {TrkBegDir = (swx*swy-sw*swyx)/(swx*swx-sw*swx2);}

                      

    // Look to add the hits
    while(CurrentPlane>0 && PlanesSinceAdded<10 && HitsAdded<4) {
                
      nhits=AllHitBank[CurrentPlane].size();
                
      if(nhits>0) {
        CurrentZPos=AllHitBank[CurrentPlane][0]->GetZPos();

        HitCam* BestHit = 0;
        MinDistanceToHit=3.*StripWidth;

        for(unsigned int j=0; j<nhits; ++j) {
          // Don't add a hit that is already in another track in the slice
          ConsiderHit=true;
          for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
            if(AllHitBank[CurrentPlane][j]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
          }
          if(ConsiderHit==false) {continue;}


          // Compare transverse positions
          TargetTPos=TrkBegTPos+TrkBegDir*(CurrentZPos-TrkBegZPos);

          centre=AllHitBank[CurrentPlane][j]->GetTPos();
          extrem1=centre+(0.055*TrkBegDir);
          extrem2=centre-(0.055*TrkBegDir);
                    
          DistanceToCurrentHit=min(fabs(centre-TargetTPos),
                                   min(fabs(extrem1-TargetTPos),fabs(extrem2-TargetTPos)));

          if(DistanceToCurrentHit<MinDistanceToHit) {
            MinDistanceToHit=DistanceToCurrentHit; 
            BestHit=AllHitBank[CurrentPlane][j]; 
            PlanesSinceAdded=0;
          }
        }

        if(BestHit) {Trk->AddHit(BestHit);}
      }
      else {PlanesSinceAdded+=2;}
      CurrentPlane-=2;
    }
  }




  // If we want to look forward into SM 2
  if(Trk->GetEndPlane()>=242 && Trk->GetEndPlane()<250) {
        
    if(PlaneView==0) {CurrentPlane=251;} else if(PlaneView==1) {CurrentPlane=250;}
    PlanesSinceAdded=0;
    TrkEndZPos=Trk->GetEndZPos();
    TrkEndTPos=Trk->GetEndTPos();
        
    
    // Get end direction
    z=0.; t=0.; sw=0.; swx=0.; swx2=0.; swy=0.; swyx=0.;
     
    // Include last few hits in track
    for(unsigned int i=0; i<Trk->GetEntries(); ++i) {
      HitCam* hit = Trk->GetHit(i);

      if(hit->GetPlane()>Trk->GetEndPlane()-8) {
        z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();
        sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z; 
      }
    }
    // Use result so far to define a tolerance
    if((swx*swx-sw*swx2)!=0) {TrkEndDir = (swx*swy-sw*swyx)/(swx*swx-sw*swx2); Tolerance = 0.2 + fabs(TrkEndDir);}
              

    // Look at some of the relevant hits across the gap
    for(int i=CurrentPlane; i<CurrentPlane+8; i+=2) {
      for(unsigned int j=0; j<AllHitBank[i].size(); ++j) {
        HitCam* hit = AllHitBank[i][j];
          
        z=hit->GetZPos(); t=hit->GetTPos(); w=hit->GetCharge();

        // If hit is within tolerance region, include with weight 0.5
        if(fabs(t-(TrkEndTPos+TrkEndDir*(z-TrkEndZPos)))<Tolerance) {
          sw+=w; swx+=w*z; swx2+=w*z*z; swy+=w*t; swyx+=w*t*z;  
        }
      }
    }

    // Get results of total fit
    if((swx*swx-sw*swx2)!=0) {TrkEndDir = (swx*swy-sw*swyx)/(swx*swx-sw*swx2);}
  
    

    // Look at hits in region
    while(CurrentPlane<500 && PlanesSinceAdded<10 && HitsAdded<4) {

      nhits=AllHitBank[CurrentPlane].size();
        
      if(nhits>0) {
        CurrentZPos=AllHitBank[CurrentPlane][0]->GetZPos();

        HitCam* BestHit = 0;
        MinDistanceToHit=3.*StripWidth;
                  
        for(unsigned int j=0; j<nhits; ++j) {
          // Don't add a hit that is already in another track in the slice
          ConsiderHit=true;
          for(unsigned int p=0; p<HitsInTracks[PlaneView].size(); ++p) {
            if(AllHitBank[CurrentPlane][j]==HitsInTracks[PlaneView][p]) {ConsiderHit=false; break;}
          }
          if(ConsiderHit==false) {continue;}


          // Compare transverse positions
          TargetTPos=TrkEndTPos+TrkEndDir*(CurrentZPos-TrkEndZPos);

          centre=AllHitBank[CurrentPlane][j]->GetTPos();
          extrem1=centre+(0.055*TrkEndDir);
          extrem2=centre-(0.055*TrkEndDir);
                    
          DistanceToCurrentHit=min(fabs(centre-TargetTPos),
                                   min(fabs(extrem1-TargetTPos),fabs(extrem2-TargetTPos)));

          if(DistanceToCurrentHit<MinDistanceToHit) {
            MinDistanceToHit=DistanceToCurrentHit; 
            BestHit=AllHitBank[CurrentPlane][j]; 
            PlanesSinceAdded=0;
          }
        }
        
        if(BestHit) {Trk->AddHit(BestHit);}
      }
      else {PlanesSinceAdded+=2;}
      CurrentPlane+=2;
    }
  }

  return;
}



void AlgTrackCamList::ClearUp()
{
  // Reset containers and tidy memory usage
  unsigned int i, j;

  TempTrack[0].clear(); TempTrack[1].clear();
  ViewSegBank[0].clear(); ViewSegBank[1].clear();
  HitsInTracks[0].clear(); HitsInTracks[1].clear(); 

  for (i=0; i<2; ++i) {
    for(j=0; j<FinalTrackBank[i].size(); ++j) {delete FinalTrackBank[i][j];}
    FinalTrackBank[i].clear();
  }


  for (i=0; i<500; ++i) {
    for(j=0; j<AllHitBank[i].size(); ++j) {if(AllHitBank[i][j]) delete AllHitBank[i][j];}
    AllHitBank[i].clear();
    HitBank[i].clear();

    for(j=0; j<ClusterBank[i].size(); ++j) {if(ClusterBank[i][j]) delete ClusterBank[i][j];}
    ClusterBank[i].clear();

    for(j=0; j<ClusterList[i].size(); ++j) {if(ClusterList[i][j]) delete ClusterList[i][j];}
    ClusterList[i].clear();

    for(j=0; j<SegmentBank[i].size(); ++j) {if(SegmentBank[i][j]) delete SegmentBank[i][j];}
    SegmentBank[i].clear();

    for(j=0; j<NDSegmentBank[i].size(); ++j) {if(NDSegmentBank[i][j]) delete NDSegmentBank[i][j];}
    NDSegmentBank[i].clear();
  }

}



void AlgTrackCamList::MatchUV(TrackCam* trku, TrackCam* trkv)
{
  // Initial check and declarations
  if(PECut<1.) {return;}

  vector<HitCam*> UTrack[500];
  vector<HitCam*> VTrack[500];


  // Store track hits  
  for(unsigned int i=0; i<trku->GetEntries(); ++i) {
    HitCam* hit = trku->GetHit(i);
    UTrack[hit->GetPlane()].push_back(hit);
  }

  for(unsigned int i=0; i<trkv->GetEntries(); ++i) {
    HitCam* hit = trkv->GetHit(i);
    VTrack[hit->GetPlane()].push_back(hit);
  }


  // Find levels of overlap
  int MinUPlane=500, MaxUPlane=-20;
  int MinUPlaneGTPECut=500, MaxUPlaneGTPECut=-20;

  int MinVPlane=500, MaxVPlane=-20;
  int MinVPlaneGTPECut=500, MaxVPlaneGTPECut=-20;

  bool GTPECut=false, ConsiderHit=true, LowPH=false;
  double PredictedTPos=0; int LastPlaneAdded;

  for(int i=0; i<500; ++i) {
    // For U
    if(UTrack[i].size()>0) {
      if(i<MinUPlane) {MinUPlane=i;}
      if(i>MaxUPlane) {MaxUPlane=i;}

      GTPECut=false;
      
      for(unsigned int j=0; j<UTrack[i].size(); ++j) {
        if(UTrack[i][j]->GetCharge()>PECut) {GTPECut=true; break;}
      }

      if(GTPECut) {
        if(i<MinUPlaneGTPECut) {MinUPlaneGTPECut=i;}
        if(i>MaxUPlaneGTPECut) {MaxUPlaneGTPECut=i;}    
      }
    }

    // For V
    if(VTrack[i].size()>0) {
      if(i<MinVPlane) {MinVPlane=i;}
      if(i>MaxVPlane) {MaxVPlane=i;}

      GTPECut=false;

      for(unsigned int j=0; j<VTrack[i].size(); ++j) {
        if(VTrack[i][j]->GetCharge()>PECut) {GTPECut=true; break;}
      }

      if(GTPECut) {
        if(i<MinVPlaneGTPECut) {MinVPlaneGTPECut=i;}
        if(i>MaxVPlaneGTPECut) {MaxVPlaneGTPECut=i;}
      }
    }
  }




  // Carry out the PH dependent UV match
  //
  // Two cases:
  // 1. Have extended too far in one view, including low PH strips
  //
  // 2. Have missed hits in one view and need to add them

  // Beginning
  if(abs(MinUPlane-MinVPlane)>8) {

    if(MinUPlane<MinVPlane) {
      // Removing hits from beginning of u track
      if(MinVPlane-MinUPlaneGTPECut<8) {
        
        for(int i=MinUPlane; i<MinUPlaneGTPECut; ++i) {
          for(unsigned int j=0; j<UTrack[i].size(); ++j) {
            UTrack[i][j]->SetUID(-1);
          }
        }
      }


      // Adding hits to beginning of v track
      else {
        LastPlaneAdded=MinVPlane;

        for(int i=MinVPlane-1; i>=MinUPlane-6; --i) {
          if(i<0 || i>120 || LastPlaneAdded-i>30) {break;}

          if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==1) {

            PredictedTPos=trkv->GetBegTPos()+trkv->GetBegDir()*(AllHitBank[i][0]->GetZPos()-trkv->GetBegZPos());

            HitCam* hit = 0;
            if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
            else {hit=AllHitBank[i][0]; LowPH=true;}

            // Don't add a hit that is already in another track in the slice
            ConsiderHit=true;
            for(unsigned int p=0; p<HitsInTracks[1].size(); ++p) {
              if(hit==HitsInTracks[1][p]) {ConsiderHit=false; break;}
            }

            // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
            if(ConsiderHit) {
              for(int q=-2; q<=2; ++q) {
                if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
                else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
              }
            }

            if(ConsiderHit==false) {continue;}

            if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trkv->AddHit(hit); LastPlaneAdded=i;}
          }      
        }
      }
    }

    else {
      // Removing hits from beginning of V track
      if(MinUPlane-MinVPlaneGTPECut<8) {
        
        for(int i=MinVPlane; i<MinVPlaneGTPECut; ++i) {
          for(unsigned int j=0; j<VTrack[i].size(); ++j) {
            VTrack[i][j]->SetUID(-1);
          }
        }
      }


      // Adding hits to beginning of u track
      else {
        LastPlaneAdded=MinUPlane;

        for(int i=MinUPlane-1; i>=MinVPlane-6; --i) {
          if(i<0 || i>120 || LastPlaneAdded-i>30) {break;}

          if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==0) {

            PredictedTPos=trku->GetBegTPos()+trku->GetBegDir()*(AllHitBank[i][0]->GetZPos()-trku->GetBegZPos());

            HitCam* hit = 0;
            if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
            else {hit=AllHitBank[i][0]; LowPH=true;}

            // Don't add a hit that is already in another track in the slice
            ConsiderHit=true;
            for(unsigned int p=0; p<HitsInTracks[0].size(); ++p) {
              if(hit==HitsInTracks[0][p]) {ConsiderHit=false; break;}
            }

            // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
            if(ConsiderHit) {
              if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
              else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
            }

            if(ConsiderHit==false) {continue;}

            if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trku->AddHit(hit); LastPlaneAdded=i;}
          }
        }
      }
    }
  }


  // End
  if(abs(MaxUPlane-MaxVPlane)>8) {

    if(MaxUPlane>MaxVPlane) {
      // Removing hits from end of u track
      if(MaxUPlaneGTPECut-MaxVPlane<8) {
        
        for(int i=MaxUPlaneGTPECut+1; i<=MaxUPlane; ++i) {
          for(unsigned int j=0; j<UTrack[i].size(); ++j) {
            UTrack[i][j]->SetUID(-1);
          }
        }
      }
      

      // Adding hits to end of v track
      else {
        LastPlaneAdded=MaxVPlane;

        for(int i=MaxVPlane; i<MaxUPlane+6; ++i) {
          if(i<0 || i>120 || i-LastPlaneAdded>30) {break;}

          if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==1) {

            PredictedTPos=trkv->GetEndTPos()+trkv->GetEndDir()*(AllHitBank[i][0]->GetZPos()-trkv->GetEndZPos());

            HitCam* hit = 0;
            if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
            else {hit=AllHitBank[i][0]; LowPH=true;}

            // Don't add a hit that is already in another track in the slice
            ConsiderHit=true;
            for(unsigned int p=0; p<HitsInTracks[1].size(); ++p) {
              if(hit==HitsInTracks[1][p]) {ConsiderHit=false; break;}
            }

            // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
            if(ConsiderHit) {
              if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
              else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
            }

            if(ConsiderHit==false) {continue;}

            if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trkv->AddHit(hit); LastPlaneAdded=i;}
          }      
        }
      }
    }
    
    else {
      // Removing hits from end of v track
      if(MaxVPlaneGTPECut-MaxUPlane<8) {
        
        for(int i=MaxVPlaneGTPECut+1; i<=MaxVPlane; ++i) {
          for(unsigned int j=0; j<VTrack[i].size(); ++j) {
            VTrack[i][j]->SetUID(-1);
          }
        }
      }
      

      // Adding hits to end of u track
      else {
        LastPlaneAdded=MaxUPlane;

        for(int i=MaxUPlane; i<MaxVPlane+6; ++i) {
          if(i<0 || i>120 || i-LastPlaneAdded>30) {break;}

          if((AllHitBank[i].size()==1 || HitBank[i].size()==1) && AllHitBank[i][0]->GetPlaneView()==0) {

            PredictedTPos=trku->GetEndTPos()+trku->GetEndDir()*(AllHitBank[i][0]->GetZPos()-trku->GetEndZPos());

            HitCam* hit = 0;
            if(HitBank[i].size()==1) {hit=HitBank[i][0]; LowPH=false;}
            else {hit=AllHitBank[i][0]; LowPH=true;}

            // Don't add a hit that is already in another track in the slice
            ConsiderHit=true;
            for(unsigned int p=0; p<HitsInTracks[0].size(); ++p) {
              if(hit==HitsInTracks[0][p]) {ConsiderHit=false; break;}
            }

            // Check to see if it is genuinely in a clean +/-5 plane instrumented part of detector
            if(ConsiderHit) {
              if(!LowPH) {for(int q=-1; q<=1; ++q) {if(q!=0 && HitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
              else {for(int q=-1; q<=1; ++q) {if(q!=0 && AllHitBank[i+2*q].size()>0) {ConsiderHit=false; break;}}}
            }

            if(ConsiderHit==false) {continue;}

            if(fabs(hit->GetTPos()-PredictedTPos)<0.4) {trku->AddHit(hit); LastPlaneAdded=i;}
          }      
        }
      }
    }
  }


  // Clear up and return
  for(int i=0; i<500; ++i) {UTrack[i].clear(); VTrack[i].clear();}
  
  return;
}



void AlgTrackCamList::Trace(const char * /* c */) const
{
}

---