StitchAlg.cxx

Go to the documentation of this file.

//
// StitchAlg.cxx
//
// echurch@fnal.gov
//
// This alg is called by TrackStich_module.cc in order to stitch together
// tracks which point at each other and have gaps to within user-set
// tolerances.


#include "larreco/RecoAlg/StitchAlg.h"

// C/C++ standard libraries
#include <time.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <cstdlib>
#include <vector>

//Framework includes:
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"

// LArSoft libraries
#include "larcoreobj/SimpleTypesAndConstants/geo_types.h" // geo::kZ


trkf::StitchAlg::StitchAlg(fhicl::ParameterSet const& pset) 
{
  ftNo = 0;
  ftListHandle.clear();
  this->reconfigure(pset); 
}

//----------------------------------------------------------

//----------------------------------------------------------
void trkf::StitchAlg::reconfigure(fhicl::ParameterSet const& pset)
{

  fCosAngTol           = pset.get< double >("CosAngTolerance", 0.95); 
  fSepTol              = pset.get< double >("SpptSepTolerance", 10.0); //cm 
}


void trkf::StitchAlg::FindHeadsAndTails( const art::Event& EvtArg, const std::string& trackModuleLabelArg  )
{

 fTrackVec.clear();
 fTrackComposite.clear();
 fHT.clear();

  EvtArg.getByLabel(trackModuleLabelArg,ftListHandle);

    // An element of fh and ft for each outer track. Keep the cos and sep parameters of the match and a string that indicates whether it's the second track's head or tail that gives the match, along with ii, jj, the indices of the outer and inner tracks.
    ft.clear(); 
    fh.clear();

    int ntrack = ftListHandle->size();
    //    std::cout << "StitchAlg.FindHeadsAndTails: Number of tracks in " << ntrack << std::endl;
    for(int ii = 0; ii < ntrack; ++ii) {
      art::Ptr<recob::Track> ptrack1(ftListHandle, ii);
      const recob::Track& track1 = *ptrack1;
      const TVector3 start1(track1.Vertex());
      const TVector3 end1(track1.End());
      const TVector3 start1Dir(track1.VertexDirection());
      const TVector3 end1Dir(track1.EndDirection());
    // For each outer track, make a vector of 1 candidate track. Doesn't need to be a vector except for desire to have a 2-iteration history.
      std::vector< std::tuple< std::string, int, int, double, double> > headvv;
      std::vector< std::tuple< std::string, int, int, double, double> > tailvv;
      
      // For head/tail keep a vector of candidate (cos,sep)
      std::vector< std::vector<std::pair< double, double>> > matchhead;
      std::vector< std::vector<std::pair< double, double>> > matchtail;
      // Neither end of track1 is yet matched:
      bool head(false);
      bool tail(false);
      std::vector<art::Ptr <recob::Track> >::iterator ith, itt;

      for(int jj = ii+1; jj < ntrack; ++jj) {
        art::Ptr<recob::Track> ptrack2(ftListHandle, jj);
        const recob::Track& track2 = *ptrack2;
        const TVector3& start2(track2.Vertex());
        const TVector3& end2(track2.End());
        const TVector3& start2Dir(track2.VertexDirection());
        const TVector3& end2Dir(track2.EndDirection());
        std::string sHT2("NA"); // track2 (receptor track) H or T is tagged as matched


        bool c12((std::abs(start1Dir.Dot(end2Dir))>fCosAngTol) && ((start1-end2).Mag()<fSepTol));
        bool c21((std::abs(end1Dir.Dot(start2Dir))>fCosAngTol) && ((start2-end1).Mag()<fSepTol));
        bool c11((std::abs(start1Dir.Dot(start2Dir))>fCosAngTol) && ((start1-start2).Mag()<fSepTol));
        bool c22((std::abs(end1Dir.Dot(end2Dir))>fCosAngTol) &&  ((end1-end2).Mag()<fSepTol));

        if ( c12 || c21 || c11 || c22 )
          {

            sHT2 = "NA";
            if (c12||c11) { head=true; }
            if (c11) { sHT2 = "H"; } else if (c12) { sHT2 = "T"; }
            if (c21||c22) { tail=true; }
            if (c21) { sHT2 = "H"; } else if (c22) { sHT2 = "T"; }

            if (head && tail) // split the tie by distance
              { 
                head = false; tail = false;
                if ( ((start1-end2).Mag() < (start2-end1).Mag()) || ((start1-end2).Mag() < (start2-end2).Mag()) || ((start1-start2).Mag() < (start2-end1).Mag()) || ((start1-start2).Mag() < (start2-end2).Mag()) )
                  { head = true; tail = false; }
                else
                  { head = false; tail = true; }
              }

            if (head) 
              {
                // 2-deep vector, for head and tail of 2nd track
                std::vector< std::pair <double,double> > headv;
                headv.push_back(std::pair<double,double>(abs(start1Dir.Dot(start2Dir)), (start1-start2).Mag()) );
                headv.push_back(std::pair<double,double>(abs(start1Dir.Dot(end2Dir)), (start1-end2).Mag()) );

                matchhead.push_back( headv );
                // if inferior, drop the new elements; if superior, replace the old
                if ( ((matchhead.size() > 1) && ( (matchhead.back().at(0).second < matchhead.front().at(0).second) || (matchhead.back().at(1).second < matchhead.front().at(1).second ) ) ) || matchhead.size()==1 )
                  {
                    if (matchhead.size()>1) matchhead.erase(matchhead.begin());
                    if (headvv.size()>1) headvv.erase(headvv.begin());
                    if (!sHT2.compare("H"))
                      {
                        auto tupTmp =std::make_tuple(sHT2,ii,jj,matchhead.back().at(0).first,matchhead.back().at(0).second);
                        headvv.push_back (tupTmp);
                      }
                    else
                      {
                        auto tupTmp = std::make_tuple(sHT2,ii,jj,matchhead.back().at(1).first,matchhead.back().at(1).second);
                        headvv.push_back (tupTmp);
                      }
                  }
                else 
                    matchhead.pop_back();

                  

                //              std::cout << "TrackStitcher: satisfied head. " << std::endl;
              } // head
                 
            else if (tail) // else if to prevent same stitching candidate being 
                               // allowed at both head and tail
              {
                // 2-deep vector, for head and tail of 2nd track
                std::vector< std::pair <double,double> > tailv;
                tailv.push_back(std::pair<double,double>( abs(end1Dir.Dot(start2Dir)),(start2-end1).Mag() ) );
                tailv.push_back(std::pair<double,double>( abs(end1Dir.Dot(end2Dir)),(end1-end2).Mag() ) );
                matchtail.push_back( tailv );
                // if inferior, drop the new elements; if superior, replace the old
                if ( ((matchtail.size() > 1) && ( (matchtail.back().at(0).second < matchtail.front().at(0).second) || (matchtail.back().at(1).second < matchtail.front().at(1).second ) ) ) || matchtail.size()==1 )
                  {
                    if (matchtail.size()>1) matchtail.erase(matchtail.begin());
                    if (tailvv.size()>1) tailvv.erase(tailvv.begin());
                    if (!sHT2.compare("T"))
                      {
                        auto tupTmp = std::make_tuple(sHT2,ii,jj,matchtail.back().at(0).first,matchtail.back().at(0).second);
                        tailvv.push_back(tupTmp);
                      }
                    else
                      {
                        auto tupTmp = std::make_tuple(sHT2,ii,jj,matchtail.back().at(1).first,matchtail.back().at(1).second);
                        tailvv.push_back(tupTmp);
                      }
                  }
                else 
                  matchtail.pop_back();

                //              std::cout << "TrackStitcher: satisfied tail. " << std::endl;
              } //tail
            /*
            std::cout << "abs(start1Dir.Dot(end2Dir)) " << std::abs(start1Dir.Dot(end2Dir)) << ", start1-end2.Mag(): " << (start1-end2).Mag() << std::endl;
            std::cout << "abs(end1Dir.Dot(start2Dir)) " << std::abs(end1Dir.Dot(start2Dir)) << ", start2-end1.Mag(): " << (start2-end1).Mag() << std::endl;
            std::cout << "abs(start1Dir.Dot(start2Dir)) " << std::abs(start1Dir.Dot(start2Dir)) << ", start1-start2.Mag(): " << (start1-start2).Mag() << std::endl;
            std::cout << "abs(end1Dir.Dot(end2Dir)) " << std::abs(end1Dir.Dot(end2Dir)) << ", end1-end2.Mag(): " << (end1-end2).Mag() << std::endl;
            std::cout << "sHT2 " << sHT2 << std::endl;      
            */
          } // end c11||c12||c21||c22

        // We've been careful to pick the best jj match for this iith track head and tail.
        // Now we need to be sure that for the jjth track head/tail we don't have two ii trks.
        if (headvv.size())
          {
            int otrk = std::get<2>(headvv.back()); // jj'th track for this iith trk
            // H or T of this jj'th trk we're matched to.
            std::string sotrkht(std::get<0>(headvv.back()));
            for (int kk=0;kk<ii;++kk)
              {
                if (std::get<2>(fh.at(kk)) == otrk && !sotrkht.compare(std::get<0>(fh.at(kk)) ) )
                  {
                    // check matching sep and pick the best one. Either erase this
                    // headvv (and it'll get null settings later below) or null out 
                    // the parameters in fh.
                    if (std::get<4>(headvv.back()) < std::get<4>(fh.at(kk)) && std::get<4>(headvv.back())!=0.0)
                      {
                        auto tupTmp2 = std::make_tuple(std::string("NA"),kk,-12,0.0,0.0);
                        fh.at(kk) = tupTmp2;
                      }
                    else if (std::get<4>(headvv.back())!=0.0)
                      {
                        headvv.pop_back();
                        break;
                      }
                  }
              }
          }
        if (tailvv.size())
          {
            int otrk = std::get<2>(tailvv.back());  // jj'th track for this iith trk
            // H or T of this jj'th trk we're matched to.
            std::string sotrkht(std::get<0>(tailvv.back()));
            for (int kk=0;kk<ii;++kk)
              {
                if (std::get<2>(ft.at(kk)) == otrk && !sotrkht.compare(std::get<0>(ft.at(kk)) ) )
                  {
                    // check matching sep and pick the best one. erase either this
                    // tailvv or null out the parameters in ft.
                    if (std::get<4>(tailvv.back()) < std::get<4>(ft.at(kk)) && std::get<4>(tailvv.back())!=0.0)
                      {
                        auto tupTmp2 = std::make_tuple(std::string("NA"),kk,-12,0.0,0.0);
                        ft.at(kk) = tupTmp2;
                      }
                    else if (std::get<4>(tailvv.back())!=0.0)
                      {
                        tailvv.pop_back();
                        break;
                      }
                  }
              }
          }


      } // jj
      
      
      auto tupTmp2 = std::make_tuple(std::string("NA"),ii,-12,0.0,0.0);
      // We always have our best 1-element tailvv and headvv for trk o at this point
      if (!headvv.size()) headvv.push_back(tupTmp2);
      if (!tailvv.size()) tailvv.push_back(tupTmp2);
      //      std::cout << "StitchAlg::FindHeadsAndTails: headvv, tailvv .get<0> is " << std::get<0>(headvv.back()) << ", " << std::get<0>(tailvv.back()) << std::endl;
      fh.push_back(headvv.back());
      ft.push_back(tailvv.back());
    
    } // ii

    //    std::cout << "fh.size, ft.size are " << fh.size() << ", " << ft.size() << std::endl;

}

void trkf::StitchAlg::FirstStitch(const std::vector<art::PtrVector <recob::Track>>::iterator itvvArg, const std::vector <recob::Track>::iterator itvArg)
{    
    // take the vector of tracks, walk through each track's vectors of xyz, dxdydz, etc 
    // and concatenate them into longer vectors. Use those to instantiate one new 
    // Stitched-together track.
    std::vector<TVector3> xyz;
    std::vector<TVector3> dxdydz;
    std::vector<TMatrixT<double> > cov;
    std::vector<double> mom;
    std::vector< std::vector <double> > dQdx;
    //art::PtrVector<recob::Track>::const_iterator


    size_t cnt(0);
    for (auto it = (*itvvArg).begin(); it!=(*itvvArg).end(); ++it)
      {

        cnt++;
        //      std::cout << "Stitching track cnt is: " << cnt << std::endl; 

        for (size_t pt = 0; pt!=(*it).get()->NumberTrajectoryPoints(); pt++)
          {
            size_t ptHere(pt);
            // ask if 1st character of 2-character string is an H. If so, reverse order
            // of the concatenation. Note however that I've had my notion of head & tail backwards
            // throughout this whole class! start1, end1 are really T, H, not H, T as I've labeled 'em till now.
            // hence flip the direction if this character is a T/H, when expecting H/T! -- EC, 29-June-2014.

            auto itvfHT = fHT.begin() + size_t (itvArg - fTrackVec.begin());
            if ( (*itvfHT).size() &&
                 (
                  // was fHT.back()
                  (cnt==1 && !(*itvfHT).at(cnt-1).compare(0,1,"H")) ||   
                   (cnt>1  && !(*itvfHT).at(cnt-2).compare(1,1,"T"))
                  )  
                 )
              ptHere = (*it).get()->NumberTrajectoryPoints() - pt - 1;

            try 
              { 
                xyz.push_back((*it).get()->LocationAtPoint(ptHere));
                //              std::cout << "Stitching track number " << cnt << " with TrajPt at ptHere " << ptHere << " at x,y,z: " << xyz.back().X() << ", " << xyz.back().Y() << ", " << xyz.back().Z() << std::endl;
                dxdydz.push_back((*it).get()->DirectionAtPoint(ptHere));
                TMatrixT<double>  dumc(5,5); 
                if (ptHere<(*it).get()->NumberCovariance())
                  dumc = (*it).get()->CovarianceAtPoint(ptHere);
                cov.push_back(dumc);
                double dumm(0.0); 
                if ((*it).get()->HasMomentum())
                  dumm = (*it).get()->MomentumAtPoint(ptHere);
                mom.push_back(dumm);
                std::vector <double> dum; 
                if (ptHere<(*it).get()->NumberdQdx(geo::kZ))
                  dum.push_back((*it).get()->DQdxAtPoint(ptHere,geo::kZ));
                else
                  dum.push_back(0.0);
                dQdx.push_back(dum);

              }
            catch (cet::exception &e)
              {
                mf::LogVerbatim("TrackStitcher bailing. ") << " One or more of xyz, dxdydz, cov, mom, dQdx elements from original Track is out of range..." << e.what() << __LINE__;
                break;
              }
          }
      }
   

    const recob::Track t(xyz,dxdydz,cov,dQdx,mom,ftNo++);
    //const art::Ptr<recob::Track> t(xyz,dxdydz,cov,dQdx,mom,ftNo++);
    fTrackVec.insert(itvArg,t);

}

void trkf::StitchAlg::WalkStitch()
{

    art::PtrVector<recob::Track> compTrack;
    std::vector <std::string> trackStatus (fh.size(),"NotDone"); // H or T 
    std::vector <std::string> HT2 ; // H or T 


    for (unsigned int ii=0; ii<fh.size(); ++ii) // same as t.size()

      {
        if (!trackStatus.at(ii).compare("Done")) continue;

        const art::Ptr<recob::Track> th(ftListHandle, ii);
        compTrack.push_back(th);
        // Should there not be an HT.push_back here??

        // start with track 1: see if head goes anywhere, walk with it. Add to compTrack.
        // Go until the other tuple's other vtx string says "NA." Then change status string 
        // of vtxsJoined to "Done" for that track.
        bool chain(true);
        int walk(ii);
        std::string sh(std::get<0>(fh.at(walk)));
        std::string st("NA");
        while (chain)
          {
            int hInd = -12; 
            int tInd = -12; 
            if (walk!=(int)ii)
              {
                sh = std::get<0>(fh.at(walk));
                st = std::get<0>(ft.at(walk));
              }

            //      std::cout << "StichAlg::WalkStitch(): Inside head chain. walk(track), sh, st, " << walk << ", " << sh << ", "<<st <<", connected to tracks: " << std::get<2>(fh.at(walk))<<", "  << std::get<2>(ft.at(walk)) <<std::endl;
            if (!sh.compare("H") || !sh.compare("T"))
              {

                hInd = std::get<2>(fh.at(walk)); // index of track that walk is connected to.
                //              std::cout << "WalkStitch(): hInd is " << hInd << std::endl;

                const art::Ptr<recob::Track> th2( ftListHandle, hInd );
                compTrack.push_back(th2);
                HT2.push_back("H"+sh);
              }

            if ((!st.compare("H") || !st.compare("T")))
              {
                tInd = std::get<2>(ft.at(walk));
                //              std::cout << "WalkStitch(): tInd is " << tInd << std::endl;
                const art::Ptr<recob::Track> th2( ftListHandle, tInd );
                compTrack.push_back(th2);
                HT2.push_back("T"+st); // since we will eventually read from 0th element forward
              }
            if (hInd!=-12) walk = hInd;
            if (tInd!=-12) walk = tInd;
            if (!sh.compare("NA") && !st.compare("NA")) 
                chain = false;

            trackStatus.at(walk) = "Done";
          } // while

        // It is possible that our first (ii'th) track had a head _and_ a tail match. Thus, we must 
        // see if tail goes anywhere. walk with it. Insert, don't push_back, to compTrack.
        chain = true;
        walk = ii;
        sh = "NA";
        st = std::get<0>(ft.at(walk));
        while (chain)
          {
            int hInd = -12; 
            int tInd = -12; 
            if (walk!=(int)ii)
              {
                sh = std::get<0>(fh.at(walk));
                st = std::get<0>(ft.at(walk));
              }

            //      std::cout << "StichAlg::WalkStitch(): Inside tail chain. walk(track), sh, st, " << walk << ", " << sh << ", "<<st <<", connected to tracks: " << std::get<2>(fh.at(walk))<<", "  << std::get<2>(ft.at(walk)) <<std::endl;
            if (!sh.compare("H") || !sh.compare("T"))
              {

                hInd = std::get<2>(fh.at(walk)); // index of track that walk is connected to.
                //              std::cout << "WalkStitch(): hInd is " << hInd << std::endl;

                const art::Ptr<recob::Track> th2( ftListHandle, hInd );
                compTrack.insert(compTrack.begin(),th2);
                HT2.insert(HT2.begin(),sh+"H");
              }

            if ((!st.compare("H") || !st.compare("T")))
              {
                tInd = std::get<2>(ft.at(walk));
                //              std::cout << "WalkStitch(): tInd is " << tInd << std::endl;
                const art::Ptr<recob::Track> th2( ftListHandle, tInd );
                compTrack.insert(compTrack.begin(),th2);
                HT2.insert(HT2.begin(),st+"T"); // since we will eventually read from 0th element forward
              }
            if (hInd!=-12) walk = hInd;
            if (tInd!=-12) walk = tInd;
            if (!sh.compare("NA") && !st.compare("NA")) 
                chain = false;

            trackStatus.at(walk) = "Done";
          } // while


        // inside FirstStitch() push_back onto the vec<vec> of components and the vec of stitched composite.
        if (compTrack.size()) 
          {
            // protect against stitching a component twice, as when somehow the same track has been matched to a 
            // head and a tail of another track. remove the repeat appearances of that track.
            for (auto iit = compTrack.begin(); iit!=compTrack.end();++iit)
              {
                int cjit(0);
                for (auto jit = iit+1; jit!=compTrack.end();++jit)
                  {
                    //              std::cout<< " cjit is ." << cjit << std::endl;
                    if (*iit==*jit) 
                      { 
                        //                      std::cout<< " About to do cjit erase." << std::endl;
                        compTrack.erase(jit); //std::cout<< "Survived the compTrack jit erase." << std::endl;
                        HT2.erase(HT2.begin()+cjit); //std::cout<< "Survived the HT2 cjit erase." << std::endl;
                        break;
                      }
                    cjit++;
                  }
              }
            fTrackComposite.push_back(compTrack);
            fHT.push_back(HT2);
            //      std::cout << "WalkStitch:: calling FirstStitch(). fTrackComposite.size(), fHT.size() " << fTrackComposite.size()<< ", " << fHT.size()  << std::endl;
            //std::cout << "WalkStitch:: calling FirstStitch(). fTrackComposite.back().size(), fHT.back().size() " << fTrackComposite.back().size()<< ", " << fHT.back().size()  << std::endl;
            /*
            for (unsigned int ll=0; ll<fHT.back().size() ; ++ll) 
              { 
                std::cout << fHT.back().at(ll); std::cout << ", "; 
              }
                   std::cout << "." << std::endl;
            */
            // want the last vector of fTrackComposite, and will want to insert on fTrackVec, hence
            // -1 and not -1, respectively.
            FirstStitch(fTrackComposite.end()-1, fTrackVec.end()); 
          }
        compTrack.clear();
        HT2.clear();
        
            
      } // end ii loop on head/tail vector of tuples.
        
}

    // This method will join separate composite tracks in which one common component is joined by its head in one
    // and by its tail in another. This can happen if the common component has a higher track index than either of
    // the the two it is separately stitched to. e.g., ____(1) -------(4) ___________(3).
    // 
bool  trkf::StitchAlg::CommonComponentStitch()
{
  // "os" for outer scope.
  int osciit(-12), oscjit(-12);
  std::vector < art::PtrVector<recob::Track> >::iterator osiComposite, osjComposite;
  art::PtrVector < recob::Track >::iterator osiAgg, osjAgg;
      

  bool match(false);
  int ciit(0); 
  for (auto iit = fTrackComposite.begin(); iit!=fTrackComposite.end()&&!match; ++iit)
    {
      ciit++;
      int cjit(ciit);
      for (auto jit = iit+1; jit!=fTrackComposite.end()&&!match; ++jit)
        {
          cjit++;
          for (auto iiit = iit->begin(); iiit!=iit->end()&&!match; ++iiit)
            {
              for (auto jjit = jit->begin(); jjit!=jit->end()&&!match; ++jjit)
                {
                  if (*iiit == *jjit) // 2 components from 2 different composites are the same
                    {
                      // head is attached to one trk and tail to another. 
                      //                      std::cout << "StitchAlg::CommonComponentStitch: We have two aggregate tracks that have a common component and need to be further stitched. " << std::endl;
                      
                      match = true;
                      osiComposite = iit;
                      osjComposite = jit;
                      osciit = ciit;
                      oscjit = cjit;
                      osiAgg = iiit;
                      osjAgg = jjit; // yes, unneeded, but we keep it for notational clarity
                      
                    }
                }
              
            }
        }
    }
  if (!match) return match;
      
  // Proceed to stitch 'em all together, dropping the one redundant component. Then
  // erase the first occurence of the composite and the matching aggregate trk.

          
  //  std::cout << "StitchAlg::CommonComponentStitch: pre erase: " << osiComposite->size() << std::endl;
  (*osiComposite).erase(osiAgg);                // erase redundant component track
                                                // do not erase this fHT element, however

  //  std::cout << "StitchAlg::CommonComponentStitch: post erase: " << osiComposite->size() << std::endl;
  // std::cout << "StitchAlg::CommonComponentStitch: fHT.size(): " << fHT.size() << std::endl;

  // Next is a loop over all remaining components in osiComposite.
  // insert the non-redundant osciit tracks onto front (back) of osjit
  // insert the non-redundant osciit vtx links onto front (back) of fHT.begin()+oscjit-1

  std::vector< std::vector<std::string> >::iterator siit(fHT.begin()+osciit-1);
  std::vector< std::vector<std::string> >::iterator sjit(fHT.begin()+oscjit-1);
  size_t itdiff(osiComposite->end()-osiComposite->begin());
  if (osjAgg == osjComposite->begin())
    {
      
      //      std::cout << "StitchAlg::begin insert starting: " << std::endl;
      // std::cout << "StitchAlg::CommonComponentStitch: itdiff: " << itdiff << std::endl;
      //std::cout << "StitchAlg::CommonComponentStitch: osiComposite.end-begin: " <<  osiComposite->end()-osiComposite->begin()<< std::endl;
      //std::cout << "StitchAlg::CommonComponentStitch: osjComposite.end-begin: " <<  osjComposite->end()-osjComposite->begin()<< std::endl;
      (*osjComposite).insert(osjComposite->begin(),osiComposite->begin(),osiComposite->begin()+itdiff); 
      //std::cout << "StitchAlg::CommonComponentStitch: siit.end-begin: " <<  siit->end()-siit->begin()<< std::endl;
      //std::cout << "StitchAlg::CommonComponentStitch: sjit.end-begin: " <<  sjit->end()-sjit->begin()<< std::endl;
      (*sjit).insert(sjit->begin(),siit->begin(),siit->begin()+itdiff);
      //std::cout << "StitchAlg::begin insert done: " << std::endl;
    }
  else if (osjAgg == (osjComposite->end()-1))
    {
      //      std::cout << "StitchAlg::end insert starting: " << std::endl;
      (*osjComposite).insert(osjComposite->end(),osiComposite->begin(),osiComposite->begin()+itdiff); 
      (*sjit).insert(sjit->end(),siit->begin(),siit->begin()+itdiff);
      //      std::cout << "StitchAlg::end insert done: " << std::endl;
    }

  //  std::cout << "StitchAlg:: 1: " << std::endl;  
  fTrackVec.erase(fTrackVec.begin()+oscjit-1);   // erase old Stitched Track, which we'll recreate now...

  //  std::cout << "StitchAlg:: 2: " << std::endl;  
  FirstStitch(osjComposite,fTrackVec.begin()+oscjit-1); // Create new Stitched Track
  //  fTrackComposite.insert(fTrackComposite.begin()+oscjit-1-1,*osjComposite);           // erase old composite Track
  //  std::cout << "StitchAlg:: 3: " << std::endl;  
  fTrackVec.erase(fTrackVec.begin()+osciit-1);   // erase old Stitched Track
  //  std::cout << "StitchAlg:: 6: " << std::endl;  
  fTrackComposite.erase(osiComposite);           // erase old composite Track
  //  std::cout << "StitchAlg:: 4: " << std::endl;  
  fHT.erase(fHT.begin()+osciit-1);               // erase old vec of vtx links
  //  std::cout << "StitchAlg:: 5: " << std::endl;  
  
  return match;
  
}   // end of bool CommonComponentStitch()