FullLDCTrackingAlg.cpp

            covMatSi[iCov] = trkExtSi->getCovMatrix()[iCov];                
            covMat[iCov] = covMatSi[iCov];
          }

          float scaling = sqrt(covMatTPC[5]/covMatSi[5]);
          covMat[5]  = covMatTPC[5];
          covMat[3]  = scaling*covMatSi[3];
          covMat[4]  = scaling*covMatSi[4];
          covMat[8]  = scaling*covMatSi[8];
          covMat[12] = scaling*covMatSi[12];

          OutputTrack->setCovMatrix(covMat);
          TrackerHitExtendedVec tpcHitVec = trkExtTPC->getTrackerHitExtendedVec();
          TrackerHitExtendedVec siHitVec = trkExtSi->getTrackerHitExtendedVec();              

          int nTPCHits = int( tpcHitVec.size());
          int nSiHits = int( siHitVec.size());        

          float edges[2];
          edges[0] = 1.0e+20;
          edges[1] = -1.0e+20;

          // find the max and min z extents from hits
          for (int iH=0;iH<nSiHits;++iH) {
            TrackerHitExtended * hitExt = siHitVec[iH];
            OutputTrack->addTrackerHitExtended(hitExt);
            hitExt->setUsedInFit(true);
	    edm4hep::ConstTrackerHit hit = hitExt->getTrackerHit();
            float zpos = float(hit.getPosition()[2]);
            if (zpos<edges[0])
              edges[0] = zpos;
            if (zpos>edges[1])
              edges[1] = zpos;
          }       
          for (int iH=0;iH<nTPCHits;++iH) {
            TrackerHitExtended * hitExt = tpcHitVec[iH];
            OutputTrack->addTrackerHitExtended(hitExt);
            hitExt->setUsedInFit(true); 
	    edm4hep::ConstTrackerHit hit = hitExt->getTrackerHit();
            float zpos = float(hit.getPosition()[2]);
            if (zpos<edges[0])
              edges[0] = zpos;
            if (zpos>edges[1])
              edges[1] = zpos;
          }

          group->setEdges(edges);
          OutputTrack->setChi2(diffMin); // will be replaced if necessary
          OutputTrack->setNDF(int(1));   // will be replaced if necessary

          _allCombinedTracks.push_back( OutputTrack );
          allMergedTracks.push_back( OutputTrack );

        }
      }
    
    } // end of loop over TPC tracks
    
    
    // check that there are some merged tracks to process
    int nMerged = int(allMergedTracks.size());

    if (nMerged>0) {
    
      // sort all merged tracks by Chi2/NDF
      // although due to the fact that NDF for these tracks is set to the value 1 above,
      // it is really only a sort on Chi2 which was really the weighted difference in angle and omega

      Sorting(allMergedTracks);

      // loop over all merged tracks
      for (int iM=0;iM<nMerged;++iM) {

        
        TrackExtended * mergedTrack = allMergedTracks[iM];
        GroupTracks * grpTrk = mergedTrack->getGroupTracks();
        
        TrackExtendedVec trkVec = grpTrk->getTrackExtendedVec();
        TrackExtended * trkTPC = NULL;
        TrackExtended * trkSi = NULL;

        int nT = int(trkVec.size());

        // only consider tracks which have been composed of excactly 2 sub tracks
        if (nT==2) {

          trkTPC = trkVec[0];
          trkSi = trkVec[1];

          GroupTracks * groupTPC = trkTPC->getGroupTracks();
          GroupTracks * groupSi  = trkSi->getGroupTracks();

          // check that both the TPC and SI track have not already been combined with other tracks ...
          if (groupTPC == NULL && groupSi == NULL) {

            // set the grouping, meaning that these tracks will not be considered further
            trkTPC->setGroupTracks( grpTrk );
            trkSi->setGroupTracks( grpTrk );

            TrackerHitExtendedVec hitVec = mergedTrack->getTrackerHitExtendedVec();
            
            int nhits = int(hitVec.size());

            int totNdf = 2*nhits - 5;
            float totChi2 = trkTPC->getChi2() + trkSi->getChi2();

            mergedTrack->setNDF( totNdf );
            mergedTrack->setChi2( totChi2 );

            if (_debug >= 2) {
              int iopt = 2;
              // PrintOutMerging(trkTPC,trkSi,iopt);
            }
            _trkImplVec.push_back( mergedTrack );
          }
        }
      }
    }
  } // end of _forceMerging
  
  
  
  // clear buffer vector
  allMergedTracks.clear();
  
  // merging split up TPC segments
  if (_mergeTPCSegments) {

    std::vector<GroupTracks*> TPCSegments;
    TPCSegments.clear();

    int nNonAssignedTPCSeg = 0;

    // loop over all TPC Tracks
    for (int i=0;i<nTPCTrk;++i) {

      TrackExtended * trkExt = _allTPCTracks[i];
      GroupTracks * group = trkExt->getGroupTracks();

      debug() << " *****************  AddNotCombinedTracks: Check track " << trkExt << " id = " << trkExt->getTrack().id()  << endmsg;
      
      // only consider those tracks which have not yet been combined
      if (group == NULL) {

        // find the min and max z extents using the hits
        TrackerHitExtendedVec currentVec = trkExt->getTrackerHitExtendedVec();
        int nCur = int(currentVec.size());
        float zmin = 1e+20;
        float zmax = -1e+20;

        for (int iCur=0;iCur<nCur;++iCur) {
          TrackerHitExtended * curTrkHitExt = currentVec[iCur];
	  edm4hep::ConstTrackerHit curTrkHit = curTrkHitExt->getTrackerHit();
          float zpos = float(curTrkHit.getPosition()[2]);
          if (zpos < zmin)
            zmin = zpos;
          if (zpos > zmax)
            zmax = zpos;
        }

        
        nNonAssignedTPCSeg++;

        // current number of TPC segment groupings
        int nGroups = int(TPCSegments.size());

        float dPtMin = 1.0e+10;
        GroupTracks * groupToAttach = NULL;
        TrackExtended * trkToAttach = NULL;

        // loop over the current TPC segment groupings
        for (int iG=0;iG<nGroups;++iG) {

          GroupTracks * segments = TPCSegments[iG];
          TrackExtendedVec segVec = segments->getTrackExtendedVec();

          // number of segments with the candidate group
          int nTrk = int(segVec.size());
          bool consider = true;

          
          if (_forbidOverlapInZTPC==1) { // if overlap in Z of the two segments is forbidden

            // loop over all tracks in the current grouping
            for (int iTrk=0;iTrk<nTrk;++iTrk) {

              TrackExtended * trkInGroup = segVec[iTrk];

              // get the number of hits from the track
              TrackerHitExtendedVec hitInGroupVec = trkInGroup->getTrackerHitExtendedVec();
              int nHitsInGrp = int(hitInGroupVec.size());

              // loop over the hits and make sure that there is no overlap in z
              for (int iHitInGrp=0;iHitInGrp<nHitsInGrp;iHitInGrp++) {

                TrackerHitExtended * xTrkExt = hitInGroupVec[iHitInGrp];
		edm4hep::ConstTrackerHit xTrk = xTrkExt->getTrackerHit();

                float xZ = float(xTrk.getPosition()[2]);
                if (xZ>zmin&&xZ<zmax) {
                  consider = false;
                  break;
                }
              }
              if (!consider)
                break; // if the candiate track's min and max z are within that of the group
            }
          }

          if (consider) {

            // again loop over the tracks in the current group
            for (int iTrk=0;iTrk<nTrk;++iTrk) {

              TrackExtended * trkInGroup = segVec[iTrk];
              int iComp = 1;
              
              // compare the tracks ... 
              float dPt = CompareTrk(trkExt,trkInGroup,_d0CutToMergeTPC,_z0CutToMergeTPC,iComp);

              // check that this tracks has the lowest delta pt  and vetomerge (i.e. fullfill momentum cut and makes sure a large number of hits have not been lost in the merger)
              if (dPt < dPtMin && !VetoMerge(trkExt,trkInGroup)) {

                dPtMin = dPt;
                groupToAttach = segments;
                trkToAttach = trkInGroup;
                if (_debug>=3) {
                  int iopt = 5;
                  // PrintOutMerging(trkExt,trkInGroup,iopt);
                }
              }
              else {
                if (_debug>=3) {
                  int iopt = 9;
                  // PrintOutMerging(trkExt,trkInGroup,iopt);
                }
              }
            }
          }
          else {
            if (_debug >= 3) {
              int iopt = 9;
              for (int iTrk=0;iTrk<nTrk;++iTrk) {
                TrackExtended * trkInGroup = segVec[iTrk];
                int iComp = 1;
                float dPt = CompareTrk(trkExt,trkInGroup,_d0CutToMergeTPC,_z0CutToMergeTPC,iComp);              
                if (dPt >= dPtMin) {          
                  // PrintOutMerging(trkExt,trkInGroup,iopt);
                }
              }
            }
          }
        }

        // check the pt cut for merging and that a group has been found to match ..
        if (dPtMin < _dPCutToMergeTPC && groupToAttach != NULL) {
          
          // add the track to the group
          groupToAttach->addTrackExtended(trkExt);
          trkExt->setGroupTracks(groupToAttach);

          // set the min and max z extents
          float zminGroup = groupToAttach->getEdges()[0];
          float zmaxGroup = groupToAttach->getEdges()[1];

          float edges[2];
          edges[0] = zmin;
          if (zminGroup<zmin)
            edges[0] = zminGroup;
          edges[1] = zmax;
          if (zmaxGroup>zmax)
            edges[1] = zmaxGroup;
          groupToAttach->setEdges(edges);
          if (_debug>=3) {
            int iopt = 3;
            // PrintOutMerging(trkExt,trkToAttach,iopt);
          }
        } else {
          
          // create a new group of segments 
          GroupTracks * newSegment = new GroupTracks(trkExt);
          trkExt->setGroupTracks(newSegment);
          debug() << " *****************  AddNotCombinedTracks: Create new TPC Segment Group for track " << trkExt << " id = " << trkExt->getTrack().id()  << endmsg;
          
          TPCSegments.push_back(newSegment);
          float edges[2];
          edges[0] = zmin;
          edges[1] = zmax;
          newSegment->setEdges(edges);
        }
      }
    }

    // At this stage all tpc segements will have been grouped.
    
    // Now try to combine the groups of TPC segments with the
    // reconstructed tracks which have already been combined into full tracks
    // containing both Si and TPC segments
    
    int nCombTrk = int(_trkImplVec.size());
    int nSegments = int(TPCSegments.size());

    //    std::cout << "Combined tracks = " << nCombTrk << endmsg;
    //    std::cout << "nSegments = " << nSegments << endmsg;

    // loop over all the TPC segment collections
    for (int iS=0;iS<nSegments;++iS) {

      GroupTracks * segments = TPCSegments[iS];
      TrackExtendedVec segVec = segments->getTrackExtendedVec();

      float zminTPCSeg = segments->getEdges()[0];
      float zmaxTPCSeg = segments->getEdges()[1];

      int nTrk = int(segVec.size());
      TrackExtended * CombTrkToAttach = NULL;
      TrackExtended * keyTrack = NULL;

      float deltaPtMin = _dPCutToMergeTPC;

      // search over the combined (good) tracks
      for (int iCTrk=0;iCTrk<nCombTrk;++iCTrk) {

        TrackExtended * combTrk = _trkImplVec[iCTrk];
        GroupTracks * groupComb = combTrk->getGroupTracks();

        bool consider = true;

        if (_forbidOverlapInZComb==1) { // if overlap in Z of the two segments is forbidden
          float zminComb = groupComb->getEdges()[0];
          float zmaxComb = groupComb->getEdges()[1];
          consider = (zminTPCSeg>zmaxComb) || (zmaxTPCSeg<zminComb);
        }
        
        // if there are not overlaps in z, if _forbidOverlapInZComb is set above
        if (consider) {

          // loop over the TPC segments in the group
          for (int iTrk=0;iTrk<nTrk;++iTrk) {

            TrackExtended * trk = segVec[iTrk];
            int iopt = 0;

            // test for compatibility
            float dPt = CompareTrk(trk,combTrk,_d0CutToMergeTPC,_z0CutToMergeTPC,iopt);
            float angleSignificance(0.);
            float significance = CompareTrkIII(trk,combTrk,_d0CutToMergeTPC,_z0CutToMergeTPC,iopt,angleSignificance);
 
            // check if this is a better match than any before
            if ( (dPt<deltaPtMin || significance <5 ) ) {
              if(VetoMerge(trk,combTrk)==false){

                // asign the track to be attached
                CombTrkToAttach = combTrk;
                keyTrack = trk;
                deltaPtMin = dPt;
              }
            }
            else {
              if (_debug>=3) {
                GroupTracks * groupCur = combTrk->getGroupTracks();
                TrackExtended * dummySi = groupCur->getTrackExtendedVec()[0];
                int iopt_temp = 8;
                // PrintOutMerging(trk,dummySi,iopt_temp);
              }
            }
          }
        }
        else {
          if (_debug>=3) {
            for (int iTrk=0;iTrk<nTrk;++iTrk) {
              TrackExtended * trk = segVec[iTrk];
              int iopt = 0;
              float dPt = CompareTrk(trk,combTrk,_d0CutToMergeTPC,_z0CutToMergeTPC,iopt);
              if (dPt>deltaPtMin) {
                GroupTracks * groupCur = combTrk->getGroupTracks();
                TrackExtended * dummySi = groupCur->getTrackExtendedVec()[0];
                int iopt_temp = 8;
                // PrintOutMerging(trk,dummySi,iopt_temp);
              }
            }
          }
        }
      }
      
      if (CombTrkToAttach != NULL) { // attach TPC segment to existing Comb Track
        GroupTracks * groupToAttach = CombTrkToAttach->getGroupTracks();          
        TrackExtended * SiCombTrk = groupToAttach->getTrackExtendedVec()[0];
        TrackExtended * TpcCombTrk = groupToAttach->getTrackExtendedVec()[1];

        if (_debug>=3) {
          int iopt = 4;
          // PrintOutMerging(keyTrack,SiCombTrk,iopt);
          iopt = 5;
          // PrintOutMerging(keyTrack,TpcCombTrk,iopt);      
        }

        for (int iTrk=0;iTrk<nTrk;iTrk++) {

          TrackExtended * segmentTrack = segVec[iTrk];
          groupToAttach->addTrackExtended( segmentTrack );
          segmentTrack->setGroupTracks( groupToAttach );

          TrackerHitExtendedVec hitVec = segmentTrack->getTrackerHitExtendedVec();

          int nHitSeg = int(hitVec.size());

          for (int iHS=0;iHS<nHitSeg;++iHS) {

            // take the hit from the segment and attach it to CombTrkToAttach
            // flagging it as not to be used in the fit
            TrackerHitExtended * hitExt = hitVec[iHS];
            hitExt->setUsedInFit(false);
            CombTrkToAttach->addTrackerHitExtended( hitExt );

          }
        }
      }
      else {
        if (nTrk==1) { // create a new group 
          GroupTracks * newGrp = new GroupTracks();
          segVec[0]->setGroupTracks(newGrp);
          newGrp->addTrackExtended(segVec[0]);
          TrackerHitExtendedVec TpcHitVec = segVec[0]->getTrackerHitExtendedVec();
          int nTpcH = int(TpcHitVec.size());
          for (int iTpcH=0;iTpcH<nTpcH;++iTpcH) {
            TpcHitVec[iTpcH]->setUsedInFit( true );
          }
          _trkImplVec.push_back(segVec[0]);
          _allNonCombinedTPCTracks.push_back(segVec[0]);
        }
        else { // several segments
          
          float zMin = 1.0e+20;
          TrackExtended * chosenTrack = NULL;

          // loop over the segments
          for (int iTrk=0;iTrk<nTrk;++iTrk) {

            TrackExtended * segment = segVec[iTrk];

            // get the lcio track which is behind this segemnt
	    edm4hep::ConstTrack track = segment->getTrack();
            ConstTrackerHitVec hitVec(track.trackerHits_begin(), track.trackerHits_end());

            debug() << "Group of orphaned TPC tracks: trying track " << track.id() << endmsg;
            
            int nHits = int(hitVec.size());

            // loop over it's hits
            for (int iH=0;iH<nHits;++iH) {
              
              float zPosi = fabs(hitVec[iH].getPosition()[2]);

              // if this segment has the hit closest to the IP so far
              if (zPosi<zMin) {
                // take this as the chosen track and break
                chosenTrack = segment;
                zMin = zPosi;
                break;
              }
            }
          }
          
          if (chosenTrack!=NULL) { // can't really ever be null.
	    debug() << "Group of orphaned TPC tracks: chosen track taken as " << chosenTrack->getTrack().id() << endmsg;
            
            // create a new group of tracks
            GroupTracks * newGroup = new GroupTracks();

            // first add the chosen track
            chosenTrack->setGroupTracks( newGroup );
            newGroup->addTrackExtended( chosenTrack );
            
            // loop over the segments ...
            for (int iTrk=0;iTrk<nTrk;++iTrk) {

              TrackExtended * segment = segVec[iTrk];
              
              TrackerHitExtendedVec hitVecS = segment->getTrackerHitExtendedVec();
              int nHitS = int(hitVecS.size());                  

              // loop over the hits for the current segment
              for (int iH=0;iH<nHitS;++iH) {
                TrackerHitExtended * trkHitExt = hitVecS[iH];

                if (segment!=chosenTrack) { // ... then don't add the hits to the fit
                  // set the relation between group and track
                  segment->setGroupTracks( newGroup );
                  newGroup->addTrackExtended( segment );
                  trkHitExt->setUsedInFit( false );
                  chosenTrack->addTrackerHitExtended( trkHitExt );                              
                }
                else {
                  trkHitExt->setUsedInFit( true );
                }
              }
            }
            _allNonCombinedTPCTracks.push_back(chosenTrack);
            _trkImplVec.push_back(chosenTrack);
          }
        }
      }
    }
    for (int iS=0;iS<nSegments;++iS) {
      GroupTracks * segments = TPCSegments[iS];
      delete segments;
    }
    TPCSegments.clear();
  }
  else { // adding all TPC segments to the list of tracks (track splitting is allowed)
    for (int i=0;i<nTPCTrk;++i) {
      TrackExtended * trkExt = _allTPCTracks[i];
      edm4hep::ConstTrack track = trkExt->getTrack();
      GroupTracks * group = trkExt->getGroupTracks();

      if (group == NULL) {
        debug() << " *****************  AddNotCombinedTracks: _mergeTPCSegments = " << _mergeTPCSegments << " : Add non combined TPC track " << trkExt << " id = " << track.id()  << endmsg;

        TrackerHitExtendedVec hitVec = trkExt->getTrackerHitExtendedVec();
        int nHTPC = int(hitVec.size());

        for (int iHTPC=0;iHTPC<nHTPC;++iHTPC) {
          hitVec[iHTPC]->setUsedInFit(true);
        }
        _trkImplVec.push_back(trkExt);
        _allNonCombinedTPCTracks.push_back( trkExt );

        GroupTracks * newGrp = new GroupTracks();
        newGrp->addTrackExtended( trkExt );
        trkExt->setGroupTracks( newGrp );

      }
    }    
  }
  
  for (int i=0;i<nSiTrk;++i) { // adding left-over Si segments to the list of tracks
    TrackExtended * trkExt = _allSiTracks[i];
    edm4hep::ConstTrack track = trkExt->getTrack();
    GroupTracks * group = trkExt->getGroupTracks();

    if (group == NULL) {
      debug() << " *****************  AddNotCombinedTracks: Add non combined Silicon Track : " << trkExt << " id = " << track.id()  << endmsg;
      
      TrackerHitExtendedVec hitVec = trkExt->getTrackerHitExtendedVec();
      int nHSi = int(hitVec.size());

      for (int iHSi=0;iHSi<nHSi;++iHSi) {
        hitVec[iHSi]->setUsedInFit(true);
      }

      _trkImplVec.push_back(trkExt);
      _allNonCombinedSiTracks.push_back( trkExt );

      GroupTracks * newGrp = new GroupTracks();
      newGrp->addTrackExtended( trkExt );
      trkExt->setGroupTracks( newGrp );   

    }
  }
  
}

void FullLDCTrackingAlg::CheckTracks() {  
  
  for(unsigned int i = 0; i< _trkImplVec.size();i++){
    TrackExtended *first = _trkImplVec[i];
    if(first==NULL)continue;
    float d0First = first->getD0();
    float z0First = first->getZ0();
    float omegaFirst = first->getOmega();
    float tanLFirst = first->getTanLambda();
    float phiFirst = first->getPhi();
    HelixClass helixFirst;
    helixFirst.Initialize_Canonical(phiFirst,d0First,z0First,omegaFirst,tanLFirst,_bField);
    float momFirst[3];
    momFirst[0]= helixFirst.getMomentum()[0];
    momFirst[1]= helixFirst.getMomentum()[1];
    momFirst[2]= helixFirst.getMomentum()[2];
    float pFirst    = sqrt(momFirst[0]*momFirst[0]+momFirst[1]*momFirst[1]+momFirst[2]*momFirst[2]);
    if(std::isnan(pFirst))continue;
    TrackerHitExtendedVec firstHitVec  = first->getTrackerHitExtendedVec();
    if(firstHitVec.size()<1)continue;
    
    for(unsigned int j = i+1; j<_trkImplVec.size();j++){
      TrackExtended *second = _trkImplVec[j];
      if(second==NULL)continue;
      float d0Second = second->getD0();
      float z0Second = second->getZ0();
      float omegaSecond = second->getOmega();
      float tanLSecond = second->getTanLambda();
      float phiSecond = second->getPhi();
      HelixClass helixSecond;
      helixSecond.Initialize_Canonical(phiSecond,d0Second,z0Second,omegaSecond,tanLSecond,_bField);
      float momSecond[3];
      momSecond[0] = helixSecond.getMomentum()[0];
      momSecond[1] = helixSecond.getMomentum()[1];
      momSecond[2] = helixSecond.getMomentum()[2];
      float pSecond    = sqrt(momSecond[0]*momSecond[0]+momSecond[1]*momSecond[1]+momSecond[2]*momSecond[2]);
      if(std::isnan(pSecond))continue;
      TrackerHitExtendedVec secondHitVec  = second->getTrackerHitExtendedVec();
      if(secondHitVec.size()<1)continue;
      if(firstHitVec.size()+secondHitVec.size()<10)continue;
      
      
      float pdot = (momFirst[0]*momSecond[0]+momFirst[1]*momSecond[1]+momFirst[2]*momSecond[2])/pFirst/pSecond;
      if(pdot<0.999)continue;
      // const float sigmaPOverPFirst  = sqrt(first->getCovMatrix()[5])/fabs(omegaFirst);
      // const float sigmaPOverPSecond = sqrt(second->getCovMatrix()[5])/fabs(omegaSecond);
      // const float deltaP = fabs(pFirst-pSecond);
      // const float sigmaDeltaP = sqrt(pFirst*sigmaPOverPFirst*pFirst*sigmaPOverPFirst+pSecond*sigmaPOverPSecond*pSecond*sigmaPOverPSecond);
      //      const float significance = deltaP/sigmaDeltaP;
      
      TrackExtended * combinedTrack = CombineTracks(first,second, _maxAllowedPercentageOfOutliersForTrackCombination, true);
      if(combinedTrack != NULL){
        const int minHits = std::min(firstHitVec.size(),secondHitVec.size());
        const int maxHits = std::max(firstHitVec.size(),secondHitVec.size());
        
        if( combinedTrack->getNDF() <= 2*maxHits+minHits-5){
          delete combinedTrack->getGroupTracks();
          delete combinedTrack;
          continue;
        }
        
        float d0    = combinedTrack->getD0();
        float z0    = combinedTrack->getZ0();
        float omega = combinedTrack->getOmega();
        float tanL  = combinedTrack->getTanLambda();
        float phi   = combinedTrack->getPhi();
        
        HelixClass helix;
        helix.Initialize_Canonical(phi,d0,z0,omega,tanL,_bField);
        float mom[3];
        mom[0]  = helix.getMomentum()[0];
        mom[1]  = helix.getMomentum()[1];
        mom[2]  = helix.getMomentum()[2];
        // float p = sqrt(mom[0]*mom[0]+mom[1]*mom[1]+mom[2]*mom[2]);
        // float chi2Sig =  (combinedTrack->getChi2() - combinedTrack->getNDF());
        // chi2Sig = chi2Sig/sqrt(combinedTrack->getNDF()*2);
        
        
        
        int nTpcFirst(0);
        int nUsedFirst(0);
        for(unsigned int ihit = 0;ihit<firstHitVec.size();ihit++){
          
          if( getDetectorID(firstHitVec[ihit]->getTrackerHit()) == lcio::ILDDetID::TPC) nTpcFirst++;
          
          if( firstHitVec[ihit]->getUsedInFit()==true ) nUsedFirst++;
        }
        
        
        int nTpcSecond(0);
        int nUsedSecond(0);
        for(unsigned int ihit = 0;ihit<secondHitVec.size();ihit++){
          if( getDetectorID(secondHitVec[ihit]->getTrackerHit()) == lcio::ILDDetID::TPC) ++nTpcSecond;
          if( secondHitVec[ihit]->getUsedInFit()==true ) ++nUsedSecond;
        }
        delete combinedTrack->getGroupTracks();
        delete combinedTrack;
      }
    }
  }
  
  
}


/*
 
 compare the following:
 
 i)   delta omega < 2 * _dOmegaForMerging
 ii)  delta d0    < d0Cut
 iii) delta z0    < z0Cut

 if above three cuts are passed return:
 
 the angle between the two tracks (return by reference)
   and
 the difference in omega
 
 */


float FullLDCTrackingAlg::CompareTrkII(TrackExtended * first, TrackExtended * second, 
                                              float d0Cut, float z0Cut,int iopt,float & Angle) {
  
  
  float result = 1.0e+20;
  Angle  = 1.0e+20; 
  float omegaFirst = first->getOmega();
  float omegaSecond = second->getOmega();
  float deltaOmega = fabs((omegaFirst-omegaSecond)/omegaSecond);
  if(deltaOmega> 2*_dOmegaForMerging)return result;
  
  
  float d0First = first->getD0();
  float z0First = first->getZ0();
  float d0Second = second->getD0();
  float z0Second = second->getZ0();
  
  bool isCloseInIP = (fabs(d0First-d0Second)<d0Cut);
  isCloseInIP = isCloseInIP && (fabs(z0Second-z0First)<z0Cut);
  
  
  if(!isCloseInIP)return result;
  
  float tanLFirst = first->getTanLambda();
  float phiFirst = first->getPhi();
  float tanLSecond = second->getTanLambda();
  float phiSecond = second->getPhi();
  float qFirst = PIOVER2 - atan(tanLFirst);
  float qSecond = PIOVER2 - atan(tanLSecond);
  
  Angle = (cos(phiFirst)*cos(phiSecond)+sin(phiFirst)*sin(phiSecond))*
  sin(qFirst)*sin(qSecond)+cos(qFirst)*cos(qSecond);
  Angle = acos(Angle);
  
  result = deltaOmega;
  
  return result;
  
}


/*
 
 compare the following:

 i)   significance d0    < d0Cut
 ii)  significance z0    < z0Cut
 iii) pdot > 0.999

 if above three cuts are passed return:
 
 the significance dAngle (return by reference)
 and
 the significance dOmega
 
 */

float FullLDCTrackingAlg::CompareTrkIII(TrackExtended * first, TrackExtended * second, 
                                               float d0Cut, float z0Cut,int iopt, float & AngleSignificance) {
  
  
  float result = 1.0e+20;
  
  float d0First = first->getD0();
  float z0First = first->getZ0();
  float omegaFirst = first->getOmega();
  float tanLFirst = first->getTanLambda();
  float phiFirst = first->getPhi();
  float qFirst = PIOVER2 - atan(tanLFirst);
  
  
  float d0Second = second->getD0();
  float z0Second = second->getZ0();
  float omegaSecond = second->getOmega();
  float tanLSecond = second->getTanLambda();
  float phiSecond = second->getPhi();
  float qSecond = PIOVER2 - atan(tanLSecond);
  
  
  //MB 2010 03
  float d0ErrFirst = sqrt(first->getCovMatrix()[0]);
  float z0ErrFirst = sqrt(first->getCovMatrix()[9]);
  //  float omegaErrFirst = sqrt(first->getCovMatrix()[5]);
  float phiErrFirst = sqrt(first->getCovMatrix()[2]);
  float qErrFirst = sqrt(cos(qFirst)*cos(qFirst)*first->getCovMatrix()[14]);
  //MB END
  //MB 2010 03
  float d0ErrSecond = sqrt(second->getCovMatrix()[0]);
  float z0ErrSecond = sqrt(second->getCovMatrix()[9]);
  //  float omegaErrSecond = sqrt(second->getCovMatrix()[5]);
  float phiErrSecond = sqrt(second->getCovMatrix()[2]);
  float qErrSecond = sqrt(cos(qSecond)*cos(qSecond)*second->getCovMatrix()[14]);
  //MB END
  
  
  //  bool isCloseInIP = (fabs(d0First-d0Second)<d0Cut);
  //isCloseInIP = isCloseInIP && (fabs(z0Second-z0First)<z0Cut);
  
  //MB 2010 03
  bool isCloseInIP = (fabs(d0First-d0Second)/sqrt(d0ErrFirst*d0ErrFirst+d0ErrSecond*d0ErrSecond)<d0Cut);
  isCloseInIP = isCloseInIP && (fabs(z0Second-z0First)/sqrt(z0ErrFirst*z0ErrFirst+z0ErrSecond*z0ErrSecond)<z0Cut);
  
  if (!isCloseInIP)return result;
  
  float Angle = (cos(phiFirst)*cos(phiSecond)+sin(phiFirst)*sin(phiSecond))*
  sin(qFirst)*sin(qSecond)+cos(qFirst)*cos(qSecond);
  
  
  HelixClass helixFirst;
  helixFirst.Initialize_Canonical(phiFirst,d0First,z0First,omegaFirst,tanLFirst,_bField);
  HelixClass helixSecond;
  helixSecond.Initialize_Canonical(phiSecond,d0Second,z0Second,omegaSecond,tanLSecond,_bField);
  
  float pFirst[3];
  float pSecond[3];
  float momFirst = 0;
  float momSecond = 0;
  
  for (int iC=0;iC<3;++iC) {
    pFirst[iC] = helixFirst.getMomentum()[iC];
    pSecond[iC] = helixSecond.getMomentum()[iC];
    momFirst += pFirst[iC]* pFirst[iC];
    momSecond += pSecond[iC]*pSecond[iC];
  }
  momFirst = sqrt(momFirst);
  momSecond = sqrt(momSecond);
  
  
  float pdot = (pFirst[0]*pSecond[0]+pFirst[1]*pSecond[1]+pFirst[2]*pSecond[2])/momFirst/momSecond;
  if(pdot<0.999)return result;
  
  const float sigmaPOverPFirst  = sqrt(first->getCovMatrix()[5])/fabs(omegaFirst);
  const float sigmaPOverPSecond = sqrt(second->getCovMatrix()[5])/fabs(omegaSecond);
  const float deltaP = fabs(momFirst-momSecond);
  const float sigmaPFirst = momFirst*sigmaPOverPFirst;
  const float sigmaPSecond = momSecond*sigmaPOverPSecond;
  const float sigmaDeltaP = sqrt(sigmaPFirst*sigmaPFirst+sigmaPSecond*sigmaPSecond);
  const float significance = deltaP/sigmaDeltaP;
  
  //MB 2010 03
  float errorAngle =sin(phiFirst)*sin(phiFirst)*phiErrFirst*phiErrFirst*cos(phiSecond)*cos(phiSecond)+
  sin(phiSecond)*sin(phiSecond)*phiErrSecond*phiErrSecond*cos(phiFirst)*cos(phiFirst)+
  sin(qFirst)*sin(qFirst)*qErrFirst*qErrFirst*cos(qSecond)*cos(qSecond)+
  sin(qSecond)*sin(qSecond)*qErrSecond*qErrSecond*cos(qFirst)*cos(qFirst)+
  cos(phiFirst)*cos(phiFirst)*phiErrFirst*phiErrFirst*(sin(phiSecond)*sin(qFirst)*sin(qSecond))*(sin(phiSecond)*sin(qFirst)*sin(qSecond))+
  cos(phiSecond)*cos(phiSecond)*phiErrSecond*phiErrSecond*(sin(phiFirst)*sin(qFirst)*sin(qSecond))*(sin(phiFirst)*sin(qFirst)*sin(qSecond))+
  cos(qFirst)*cos(qFirst)*qErrFirst*qErrFirst*(sin(phiFirst)*sin(phiSecond)*sin(qSecond))*(sin(phiFirst)*sin(phiSecond)*sin(qSecond))+
  cos(qSecond)*cos(qSecond)*qErrSecond*qErrSecond*(sin(phiFirst)*sin(phiSecond)*sin(qFirst))*(sin(phiFirst)*sin(phiSecond)*sin(qFirst));
  
  if(Angle<1.){
    errorAngle = sqrt(1./(1.-Angle*Angle)*errorAngle);
  }else{
    errorAngle = sqrt(errorAngle);
  }
  
  if(errorAngle<1.e-6)errorAngle=1.e-6;
  
  AngleSignificance = fabs(acos(Angle)/errorAngle);
    
  return significance;
  
}

/*
 
 compare the following:
 
 i)   delta d0    < d0Cut  and optionally (d0_1 + d0_2) < d0cut
 ii)  delta z0    < z0Cut
 
 if above two cuts are passed then:
 
 if ( (ptFirst<_PtCutToMergeTPC) && (ptSecond<_PtCutToMergeTPC) ) {
 
     then check the difference in momentum 
 
 else 
 
     check for cases where PatRec splits non-looping TPC tracks
     look for two tracks where total tpc hits are not more than total number
     of pad rows and that the hits on one track are close to the helix of the
     other track.
 
     Check that the angular and momentum difference meets the cuts for either hight or low  pt
     also check if the angle is very small between the tracks and that the significance of the difference in pt is less than 10
 
 
 .... 
 
 note currently this is only used in AddNotCombinedTracks
 
 */


float FullLDCTrackingAlg::CompareTrk(TrackExtended * first, TrackExtended * second, float d0Cut, float z0Cut,int iopt) {
  
  float result = 1.0e+20;
  
  float d0First = first->getD0();
  float z0First = first->getZ0();
  float omegaFirst = first->getOmega();
  float tanLFirst = first->getTanLambda();
  float phiFirst = first->getPhi();
  
  float d0Second = second->getD0();
  float z0Second = second->getZ0();
  float omegaSecond = second->getOmega();
  float tanLSecond = second->getTanLambda();
  float phiSecond = second->getPhi();
  
  bool isCloseInIP = (fabs(d0First-d0Second)<d0Cut);
  
  if (iopt>0) isCloseInIP = isCloseInIP || (fabs(d0First+d0Second)<d0Cut);
  
  isCloseInIP = isCloseInIP && (fabs(z0Second-z0First)<z0Cut);
  
  
  HelixClass helixFirst;
  helixFirst.Initialize_Canonical(phiFirst,d0First,z0First,omegaFirst,tanLFirst,_bField);
  HelixClass helixSecond;
  helixSecond.Initialize_Canonical(phiSecond,d0Second,z0Second,omegaSecond,tanLSecond,_bField);
  
  float pFirst[3];
  float pSecond[3];
  float dPminus[3];
  float dPplus[3];
  float momFirst = 0;
  float momSecond = 0;
  float momMinus = 0;
  float momPlus = 0;
  
  if ( isCloseInIP ) {
    
    for (int iC=0;iC<3;++iC) {

      pFirst[iC] = helixFirst.getMomentum()[iC];
      pSecond[iC] = helixSecond.getMomentum()[iC];

      momFirst += pFirst[iC]* pFirst[iC];
      momSecond += pSecond[iC]*pSecond[iC];

      dPminus[iC] = pFirst[iC] - pSecond[iC];
      dPplus[iC] = pFirst[iC] + pSecond[iC];

      momMinus += dPminus[iC]*dPminus[iC];
      momPlus += dPplus[iC]*dPplus[iC];
    }

    momFirst = sqrt(momFirst);
    momSecond = sqrt(momSecond);
    
    float ptFirst = sqrt(pFirst[0]*pFirst[0]+pFirst[1]*pFirst[1]);
    float ptSecond = sqrt(pSecond[0]*pSecond[0]+pSecond[1]*pSecond[1]);
    
    // if both track's pt are lower than _PtCutToMergeTPC
    if ( (ptFirst<_PtCutToMergeTPC) && (ptSecond<_PtCutToMergeTPC) ) {
      
      momMinus = sqrt(momMinus);
      momPlus = sqrt(momPlus);

      float nom = momMinus;

      // get the smaller difference for the nominator
      if (momPlus<nom && iopt>0)
        nom = momPlus;

      float den = momFirst;

      // get the smallest momentum for the denominator
      if (momSecond<momFirst)
        den = momSecond;
      
      result = nom/den;     
      
    }
    else {
      
      
      // check for cases where PatRec splits non-looping TPC tracks 
      // look for two tracks where total tpc hits are not more than total number
      // of pad rows and that the hits on one track are close to the helix of the
      // other track
      
      float dpOverP = 2.0*fabs(momFirst-momSecond)/(momFirst+momSecond);
      const float pdot = (pFirst[0]*pSecond[0]+pFirst[1]*pSecond[1]+pFirst[2]*pSecond[2])/momFirst/momSecond;
      const float sigmaPOverPFirst  = sqrt(first->getCovMatrix()[5])/fabs(omegaFirst);
      const float sigmaPOverPSecond = sqrt(second->getCovMatrix()[5])/fabs(omegaSecond);
      const float deltaP = fabs(momFirst-momSecond);
      const float sigmaPFirst = momFirst*sigmaPOverPFirst;
      const float sigmaPSecond = momSecond*sigmaPOverPSecond;
      const float sigmaDeltaP = sqrt(sigmaPFirst*sigmaPFirst+sigmaPSecond*sigmaPSecond);
      const float significance = deltaP/sigmaDeltaP;
      
      
      //compare angle between the two vectors (cos theta) and their momentum
      if( (pdot>_cosThetaCutHighPtMerge && dpOverP<_momDiffCutHighPtMerge) 
         || 
         (pdot>_cosThetaCutSoftHighPtMerge && dpOverP<_momDiffCutSoftHighPtMerge) 
         || (pdot > 0.9999 && significance <10) 
         ){
        
        
        int nTrkGrpFirst = 0;
        int nTrkGrpSecond = 0;
        ConstTrackerHitVec hitvecFirst;
        ConstTrackerHitVec hitvecSecond;
        GroupTracks * groupFirst = first->getGroupTracks();
        GroupTracks * groupSecond = second->getGroupTracks();

        // does the first track belong to a group ...
        // if it does then get all the hits from the group
        if(groupFirst!=NULL){
          
          TrackExtendedVec tracksInGroupFirst = groupFirst->getTrackExtendedVec();
          nTrkGrpFirst = int(tracksInGroupFirst.size());