////////////////////////////////////////////////////////////////////////
// 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 )
	//////////////                           Segm    Seg    Segp     Segp2

	// 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 
	//////////////                              Segm     Seg    Segp

	// 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 )
	//////////////                            Seg   Segp      Segp2

	// 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 
	//////////////                           Seg    Segp           

	// 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<ViewSe