TCurlFinder.cxx

Go to the documentation of this file.

//
// ... Section of Include Files
//
#include <iostream>
#include <fstream>
#include <math.h>
#include <stdlib.h>
//#include "panther/panther.h"
#include "CLHEP/String/Strings.h"
#include "TrkReco/TCurlFinder.h"
#include "TrkReco/TMDCWire.h"
#include "TrkReco/TMDCWireHit.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TCircle.h"
#include "TrkReco/TTrack.h"
#include "TrkReco/TSegmentCurl.h"
//#include "TrkReco/TSvdFinder.h"
//#include "TrkReco/TSvdHit.h"
//#include MDC_H
//#include SVD_H
 
#include "MdcTables/MdcTables.h"
 
// Following 3 parameters : (0,0,0,0) is best!
// Other cases are for the development.
#define DEBUG_CURL_DUMP 0 
#define DEBUG_CURL_SEGMENT 0 
#define DEBUG_CURL_GNUPLOT 0 
#define DEBUG_CURL_MC 0 
 
#if (DEBUG_CURL_DUMP+DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
#include <set>
#include <algo.h>
//#include "tables/belletdf.h"
static int debugMcFlag = 1;
#endif
 
//
// ... Constructor and Destructor Section ...
//
TCurlFinder::TCurlFinder(void)
  : m_builder("CurlBuilder"),
    m_fitter("TCurlFinder Fitter"),
    m_debugCdcFrame(false),
    m_debugPlotFlag(0),
    m_debugFileNumber(0) {
  // *****NOTE*****
  // Do not use this!!!!!
  // Because parameters can not be set correctly.
}
 
TCurlFinder::TCurlFinder(const unsigned min_segment,
             const unsigned min_salvage,
             const double bad_distance_for_salvage,
             const double good_distance_for_salvage,
             const unsigned min_sequence,
             const unsigned min_fullwire,
             const double range_for_axial_search,
             const double range_for_stereo_search,
             const unsigned superlayer_for_stereo_search,
             const double range_for_axial_last2d_search,
             const double range_for_stereo_last2d_search,
             const double trace2d_distance,
             const double trace2d_first_distance,
             const double trace3d_distance,
             const unsigned determine_one_track,
             //
             const double selector_max_impact,
             const double selector_max_sigma,
             const double selector_strange_pz,
             const double selector_replace_dz,
             //
             const unsigned stereo_2dfind,
             const unsigned merge_exe,
             const double merge_ratio,
             const double merge_z_diff,
             const double mask_distance,
             const double ratio_used_wire,
             const double range_for_stereo1,
             const double range_for_stereo2,
             const double range_for_stereo3,
             const double range_for_stereo4,
             const double range_for_stereo5,
                         const double range_for_stereo6,
             //
             const double z_cut,
             const double z_diff_for_last_attend,
             const unsigned svd_reconstruction,
             const double min_svd_electrons,
             const unsigned on_correction,
             const unsigned output_2dtracks,
             const unsigned curl_version,
             
             //jialk
             const  double minimum_seedLength,
             const  double minimum_2DTrackLength,
             const  double minimum_3DTrackLength,
             const  double minimum_closeHitsLength,
             const  double MIN_RADIUS_OF_STRANGE_TRACK,
             const  double ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK)
  : m_builder("CurlBuilder"),
    m_fitter("TCurlFinder Fitter"),
    m_debugCdcFrame(false),
    m_debugPlotFlag(0),
    m_debugFileNumber(0) {
  StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
  if(scmgn!=StatusCode::SUCCESS) { 
    std::cout<< __FILE__<<" Unable to open Magnetic field service"<<std::endl;
  }
  //...Set Parameter
  m_param.MIN_SEGMENT = min_segment;
  m_param.MIN_SALVAGE = min_salvage;
  m_param.BAD_DISTANCE_FOR_SALVAGE  = bad_distance_for_salvage;
  m_param.GOOD_DISTANCE_FOR_SALVAGE = good_distance_for_salvage;
  m_param.MIN_SEQUENCE = min_sequence;
  m_param.MAX_FULLWIRE = min_fullwire;
  m_param.RANGE_FOR_AXIAL_SEARCH  = range_for_axial_search;
  m_param.RANGE_FOR_STEREO_SEARCH = range_for_stereo_search;
  m_param.SUPERLAYER_FOR_STEREO_SEARCH = superlayer_for_stereo_search;
  m_param.RANGE_FOR_AXIAL_LAST2D_SEARCH  = range_for_axial_last2d_search;
  m_param.RANGE_FOR_STEREO_LAST2D_SEARCH = range_for_stereo_last2d_search;
  m_param.TRACE2D_DISTANCE = trace2d_distance;
  m_param.TRACE2D_FIRST_SUPERLAYER = trace2d_first_distance;
  m_param.TRACE3D_DISTANCE = trace3d_distance;
  m_param.DETERMINE_ONE_TRACK = determine_one_track;
  //
  m_param.SELECTOR_MAX_IMPACT = selector_max_impact;
  m_param.SELECTOR_MAX_SIGMA  = selector_max_sigma;
  m_param.SELECTOR_STRANGE_PZ = selector_strange_pz;
  m_param.SELECTOR_REPLACE_DZ = selector_replace_dz;
  //
  m_param.STEREO_2DFIND = stereo_2dfind;
  m_param.MERGE_EXE    = merge_exe;
  m_param.MERGE_RATIO  = merge_ratio;
  m_param.MERGE_Z_DIFF = merge_z_diff;
  m_param.MASK_DISTANCE = mask_distance;
  m_param.RATIO_USED_WIRE = ratio_used_wire;
  m_param.RANGE_FOR_STEREO_FIRST  = range_for_stereo1;
  m_param.RANGE_FOR_STEREO_SECOND = range_for_stereo2;
  m_param.RANGE_FOR_STEREO_THIRD  = range_for_stereo3;
  m_param.RANGE_FOR_STEREO_FORTH  = range_for_stereo4;
  m_param.RANGE_FOR_STEREO_FIFTH  = range_for_stereo5;
  m_param.RANGE_FOR_STEREO_SIXTH  = range_for_stereo6;
  //
  m_param.Z_CUT = z_cut;
  m_param.Z_DIFF_FOR_LAST_ATTEND = z_diff_for_last_attend;
  m_param.SVD_RECONSTRUCTION = svd_reconstruction;
  m_param.MIN_SVD_ELECTRONS = min_svd_electrons;
  m_param.ON_CORRECTION = on_correction;
  m_param.OUTPUT_2DTRACKS = output_2dtracks;
  m_param.CURL_VERSION = curl_version;
  m_param.now();
  //jialk
  m_param.minimum_seedLength = minimum_seedLength;
  m_param.minimum_2DTrackLength = minimum_2DTrackLength;
  m_param.minimum_3DTrackLength = minimum_3DTrackLength;
  m_param.minimum_closeHitsLength = minimum_closeHitsLength;
  m_param.MIN_RADIUS_OF_STRANGE_TRACK = MIN_RADIUS_OF_STRANGE_TRACK;
  m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK = ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK;
                                                    
 
  //...Set up TBuilder...
  m_builder.setParam(m_param);
}
 
TCurlFinder::~TCurlFinder(void) {
}
 
//
// ... Version and Name Section ...
//
std::string
TCurlFinder::name(void) const {
  return std::string("Curling Track Finder");
}
 
std::string
TCurlFinder::version(void) const {
  return std::string("3.00");
}
 
//
// ... AList Sorter Section...
//
#if defined(__GNUG__)
int
sortBySequentialLength( const TSegmentCurl **a, const TSegmentCurl **b ) {
  if( (*a)->maxSeq() < (*b)->maxSeq() ){
    return 1;
  }else if( (*a)->maxSeq() == (*b)->maxSeq() ){
    return 0;
  }else{
    return -1;
  }
}
#else
extern "C" int
sortBySequentialLength( const void *av, const void *bv ) {
  const TSegmentCurl **a((const TSegmentCurl**)av);
  const TSegmentCurl **b((const TSegmentCurl**)bv);
  if( (*a)->maxSeq() < (*b)->maxSeq() ){
    return 1;
  }else if( (*a)->maxSeq() == (*b)->maxSeq() ){
    return 0;
  }else{
    return -1;
  }
}
#endif
 
#if defined(__GNUG__)
int
sortByArcLength( const TMLink **a, const TMLink **b ) {
  if( (*a)->position().x() > (*b)->position().x() ){
    return 1;
  }else if( (*a)->position().x() == (*b)->position().x() ){
    return 0;
  }else{
    return -1;
  }
}
#else
int
sortByArcLength( const void*av, const void*bv ) {
  const TMLink **a((const TMLink**)av);
  const TMLink **b((const TMLink**)bv);
  if( (*a)->position().x() > (*b)->position().x() ){
    return 1;
  }else if( (*a)->position().x() == (*b)->position().x() ){
    return 0;
  }else{
    return -1;
  }
}
#endif
//
// ... Utility Section ...
//
double
TCurlFinder::distance(const double x, const double y) const {
  return sqrt(x*x+y*y);
}
 
unsigned
TCurlFinder::offset(const unsigned layerId) const {
  // input   - layer id#(0-49)
  // output  - offset of its layer
  if( layerId ==  0 || layerId ==  2 || layerId ==  4 ||
      layerId ==  6 || layerId ==  8 || layerId ==  9 ||
      layerId == 11 || layerId == 13 || layerId == 15 ||
      layerId == 17 || layerId == 18 || layerId == 20 ||
      layerId == 22 || layerId == 24 || layerId == 26 ||
      layerId == 27 || layerId == 29 || layerId == 31 ||
      layerId == 33 || layerId == 35 || layerId == 36 ||
      layerId == 38 || layerId == 40 || layerId == 42 ||
      layerId == 44 || layerId == 45 || layerId == 47 ||
      layerId == 49 ) return 1; // off set is 0.5
  return 0; // off set is 0
}
 
unsigned
TCurlFinder::layerId(const double &R) const {
  // R is radius for MDC but is 2*radius for track
  double r = R*10.;// cm -> mm
/*  if(r < 83.0 || r > 874.)return 50;
  if(r <=  93.0)return  0; if(r <= 103.25)return 1; if(r <= 113.5)return  2;
  if(r <= 136.0)return  3; if(r <= 152.0)return  4; if(r <= 169.0)return  5;
  if(r <= 186.0)return  6; if(r <= 202.0)return  7; if(r <= 217.0)return  8;
  if(r <= 232.0)return  9; if(r <= 248.0)return 10; if(r <= 264.0)return 11;
  if(r <= 280.0)return 12; if(r <= 297.0)return 13; if(r <= 313.0)return 14;
  if(r <= 330.0)return 15; if(r <= 346.0)return 16; if(r <= 361.0)return 17;
  if(r <= 376.0)return 18; if(r <= 392.0)return 19; if(r <= 408.0)return 20;
  if(r <= 424.0)return 21; if(r <= 441.0)return 22; if(r <= 458.0)return 23;
  if(r <= 474.0)return 24; if(r <= 490.0)return 25; if(r <= 505.0)return 26;
  if(r <= 520.0)return 27; if(r <= 536.0)return 28; if(r <= 552.0)return 29;
  if(r <= 568.0)return 30; if(r <= 585.0)return 31; if(r <= 602.0)return 32;
  if(r <= 618.0)return 33; if(r <= 634.0)return 34; if(r <= 649.0)return 35;
  if(r <= 664.0)return 36; if(r <= 680.0)return 37; if(r <= 696.0)return 38;
  if(r <= 712.0)return 39; if(r <= 729.0)return 40; if(r <= 746.0)return 41;
  if(r <= 762.0)return 42; if(r <= 778.0)return 43; if(r <= 793.0)return 44;
  if(r <= 808.0)return 45; if(r <= 824.0)return 46; if(r <= 840.0)return 47;
  if(r <= 856.0)return 48; if(r <= 874.0)return 49;
*/
//Liuqg 060915
  if(r < 73.0 || r > 772.2)return 43;  
  if(r <=  85.0)return  0; if(r <= 97.0)return 1; if(r <= 109.0)return  2;
  if(r <= 121.0)return  3; if(r <= 133.0)return  4; if(r <= 145.0)return  5;
  if(r <= 157.0)return  6; if(r <= 179.0)return  7; if(r <= 205.2)return  8;
  if(r <= 221.4)return  9; if(r <= 237.6)return 10; if(r <= 253.8)return 11;
  if(r <= 270.0)return 12; if(r <= 286.2)return 13; if(r <= 302.4)return 14;
  if(r <= 318.6)return 15; if(r <= 334.8)return 16; if(r <= 351.0)return 17;
  if(r <= 367.2)return 18; if(r <= 387.45)return 19; if(r <= 407.7)return 20;
  if(r <= 423.9)return 21; if(r <= 440.1)return 22; if(r <= 456.3)return 23;
  if(r <= 472.5)return 24; if(r <= 488.7)return 25; if(r <= 504.9)return 26;
  if(r <= 521.1)return 27; if(r <= 537.3)return 28; if(r <= 554.5)return 29;
  if(r <= 569.7)return 30; if(r <= 585.9)return 31; if(r <= 602.1)return 32;
  if(r <= 618.3)return 33; if(r <= 634.5)return 34; if(r <= 654.75)return 35;
  if(r <= 675.0)return 36; if(r <= 691.2)return 37; if(r <= 707.4)return 38;
  if(r <= 723.6)return 39; if(r <= 739.8)return 40; if(r <= 756.0)return 41;
  if(r <= 772.2)return 42;
}
 
unsigned 
TCurlFinder::maxLocalLayerId(const unsigned superLayerId) const {
  // input   - superlayer id#(0-10)
  // output  - max# id of locallayer in its superlayer.
//Liuqg 060926, change to BES, superLayerId means the max num of layers in each superlayers.
/*  if(superLayerId == 1 || superLayerId == 3)return 2;
  if(superLayerId == 5 || superLayerId == 7 ||
     superLayerId == 9)return 3;
  if(superLayerId == 4 || superLayerId == 6 ||
     superLayerId == 8 || superLayerId == 10)return 4;
  if(superLayerId == 0 || superLayerId ==  2)return 5; */
  if(superLayerId == 10)return 2;
  if(superLayerId >= 0 && superLayerId < 10)return 3;
 
  std::cerr << "Error in the CurlFinder(maxLocalLayerId). superLayerId = " 
       << superLayerId  << std::endl;
  return 0;
}
 
int
TCurlFinder::nextSuperAxialLayerId(const unsigned superLayerId, const int in) const {
  // This function is to find next axial superlayer! 
  // input   - superLayerId = superlayer id#(0-10)
  //         - in = find inner axial superlayer from "superLayerID" 
  //                This is depth of it. 
  //                ex1. If superLayerID = 2, in =  1, return 0 
  //                ex2. If superLayerID = 2, in = -1, return 4 
  // output  - return 0, 2, 4, 6, 8, 10 = no error
  //                                 -1 = error
  if(superLayerId > 10 || superLayerId < 0){
    std::cerr << "Error in the CurlFinder(nextSuperAxialLayerId)." << std::endl;
    return -1;
  }
//Liuqg 060920, the following numbers have been changed from belle to besiii  
  if(in == 0){
    if(superLayerId == 2 || superLayerId ==  3 ||
       superLayerId == 4 || superLayerId ==  9 ||
       superLayerId == 10){
      return superLayerId;
    }else{
      //return superLayerId - 1;
      return -1;
    }
  }
  // almost case --> inner type
  if((superLayerId == 3 ) && in == 1)return 2;
  if(superLayerId == 4){
    if(in == 1)return 3; if(in == 2)return 2;
  }
  if(superLayerId == 5 || superLayerId ==  6 ||
     superLayerId == 7 || superLayerId ==  8 ||
     superLayerId == 9){
    if(in == 1)return 4; if(in == 2)return 3;
    if(in == 3)return 2;
  }
  if(superLayerId == 10){
    if(in == 1)return 9; if(in == 2)return 4;
    if(in == 3)return 3; if(in == 4)return 2;
  }
  // rare case --> outer type
  if(superLayerId == 0 || superLayerId ==  1){
    if(in == -1)return  2; if(in == -2)return 3;
    if(in == -3)return  4; if(in == -4)return 9;
    if(in == -5)return  10;
  }
  if( superLayerId == 2){
    if(in == -1)return  3; if(in == -2)return 4;
    if(in == -3)return  9; if(in == -4)return 10;
  }
  if(superLayerId == 3){
    if(in == -1)return 4; if(in == -2)return 9;
    if(in == -3)return 10;
  }
  if(superLayerId == 4){
    if(in == -1)return  9; if(in == -2)return 10;
  }
  if(superLayerId == 5 || superLayerId ==  6 ||
     superLayerId == 7 || superLayerId ==  8 ||
     superLayerId == 9){
    if(in == -1)return 10;
  }
 
  return -1;
}
 
int
TCurlFinder::nextSuperStereoLayerId(const unsigned superLayerId, const int in) const {
  // This function is to find next stereo superlayer! 
  // input   - superLayerId = superlayer id#(0-10)
  //         - in = find inner stereo superlayer from "superLayerID" 
  //                This is depth of it. 
  //                ex1. If superLayerID = 2, in =  1, return 1 
  //                ex2. If superLayerID = 2, in = -1, return 3 
  // output  - return 1, 3, 5, 7, 9 = no error
  //                             -1 = error
  if(superLayerId > 10 || superLayerId < 0){
    std::cerr << "Error in the CurlFinder(nextSuperStereoLayerId)." << std::endl;
    return -1;
  }
//Liuqg 060920, the following numbers have been changed from belle to besiii  
  if(in == 0){
    if(superLayerId == 0 || superLayerId ==  1 ||
       superLayerId == 5 || superLayerId ==  6 ||
       superLayerId == 7 || superLayerId ==  8){
      return superLayerId;
    }else{
      return -1;
    }
  }
  // almost case --> inner type
  if((superLayerId == 1 || superLayerId ==  2 ||
      superLayerId == 3 || superLayerId ==  4) && in == 1)return 0;
  if(superLayerId == 5){
    if(in == 1)return 1; if(in == 2)return 0;
  }
  if(superLayerId == 6){
    if(in == 1)return 5; if(in == 2)return 1;
    if(in == 3)return 0;
  }
  if(superLayerId == 7){
    if(in == 1)return 6; if(in == 2)return 5;
    if(in == 3)return 1; if(in == 4)return 0;
  }
  if(superLayerId ==  8){
    if(in == 1)return 7; if(in == 2)return 6;
    if(in == 3)return 5; if(in == 4)return 1;
    if(in == 5)return 0;
  }
  if(superLayerId ==  9 || superLayerId ==  10){
    if(in == 1)return 8; if(in == 2)return 7;
    if(in == 3)return 6; if(in == 4)return 5;
    if(in == 5)return 1; if(in == 6)return 0;
  }
  // rare case --> outer type
  if(superLayerId == 0){
    if(in == -1)return 1; if(in == -2)return 5;
    if(in == -3)return 6; if(in == -4)return 7;
    if(in == -5)return 8;
  }
  if(superLayerId == 1 || superLayerId == 2 ||
     superLayerId == 3 || superLayerId == 4){
    if(in == -1)return 5; if(in == -2)return 6;
    if(in == -3)return 7; if(in == -4)return 8;
  }
  if(superLayerId == 5){
    if(in == -1)return 6; if(in == -2)return 7;
    if(in == -3)return 8;
  }
  if(superLayerId == 6){
    if(in == -1)return 7; if(in == -2)return 8;
  }
  if(superLayerId == 7){
    if(in == -1)return 8;
  }
 
  return -1;
}
 
unsigned 
TCurlFinder::nAxialHits(const double &r) const {
  // r is radius for MDC but is 2*radius for track
  const double eps = 0.2;
//changed to BESIII, Liuqg 060921
/*  if(r < 8.8-eps)        return 0;
  else if(r < 9.8-eps)   return 1;
  else if(r < 10.85-eps) return 2;
  else if(r < 12.8-eps)  return 3;
  else if(r < 14.4-eps)  return 4;
  else if(r < 15.95-eps) return 5;
  else if(r < 22.45-eps) return 6;
  else if(r < 24.0-eps)  return 7;
  else if(r < 25.6-eps)  return 8;
  else if(r < 27.2-eps)  return 9;
  else if(r < 28.8-eps)  return 10;
  else if(r < 30.4-eps)  return 11;
  else if(r < 36.85-eps) return 12;
  else if(r < 38.4-eps)  return 13;
  else if(r < 40.0-eps)  return 14;
  else if(r < 41.6-eps)  return 15;
  else if(r < 43.15-eps) return 16;
  else if(r < 51.25-eps) return 17;
  else if(r < 52.8-eps)  return 18;
  else if(r < 54.4-eps)  return 19;
  else if(r < 56.0-eps)  return 20;
  else if(r < 57.55-eps) return 21;
  else if(r < 65.65-eps) return 22;
  else if(r < 67.2-eps)  return 23;
  else if(r < 68.8-eps)  return 24;
  else if(r < 70.4-eps)  return 25;
  else if(r < 71.9-eps)  return 26;
  else if(r < 80.05-eps) return 27;
  else if(r < 81.6-eps)  return 28;
  else if(r < 83.2-eps)  return 29;
  else if(r < 84.8-eps)  return 30;
  else if(r < 86.3-eps)  return 31;
  else return 32;
*/
  if(r < 20.52-eps)        return 0;
  else if(r < 22.14-eps)   return 1;
  else if(r < 23.76-eps)   return 2;
  else if(r < 25.38-eps)   return 3;
  else if(r < 27.0-eps)    return 4;
  else if(r < 28.62-eps)   return 5;
  else if(r < 30.24-eps)   return 6;
  else if(r < 31.86-eps)   return 7;
  else if(r < 33.48-eps)   return 8;
  else if(r < 35.1-eps)    return 9;
  else if(r < 36.72-eps)   return 10;
  else if(r < 38.34-eps)   return 11;
  else if(r < 67.5-eps)    return 12;
  else if(r < 69.12-eps)   return 13;
  else if(r < 70.74-eps)   return 14;
  else if(r < 72.36-eps)   return 15;
  else if(r < 73.98-eps)   return 16;
  else if(r < 75.6-eps)    return 17;
  else if(r < 77.22-eps)   return 18;
}
 
//
// ... Clean or Remove Section ...
//
void 
TCurlFinder::makeList(AList<TMLink> &madeList, 
              const AList<TSegmentCurl> &originalList, 
              const AList<TMLink> &removeList) {
  // This is to make "madeList" from "originalList",
  // but remove "removeList" from this "madeList", that is,
  // madeList = originalList - removeList.
  madeList.removeAll();
  for(unsigned i = 0, size = originalList.length(); i < size; ++i)
    madeList.append(originalList[i]->list());
  madeList.remove(removeList);
}
 
void 
TCurlFinder::makeList(AList<TMLink> &madeList, 
              const AList<TMLink> &originalList, 
              const AList<TMLink> &removeList) {
  // This is to make "madeList" from "originalList",
  // but remove "removeList" from this "madeList", that is,
  // madeList = originalList - removeList.
  madeList.removeAll();
  madeList.append(originalList);
  madeList.remove(removeList);
}
 
void
TCurlFinder::clear(void) {
  // This is to clear this Class(TCurlFinder) in TrkReco.cc . 
  // Private members are cleaned.
  HepAListDeleteAll(m_allAxialHitsOriginal);
  HepAListDeleteAll(m_allStereoHitsOriginal);
  HepAListDeleteAll(m_segmentList);
  HepAListDeleteAll(m_allCircles);
  HepAListDeleteAll(m_allTracks);
  
  m_unusedAxialHitsOriginal.removeAll();
  m_unusedStereoHitsOriginal.removeAll();
  m_unusedAxialHits.removeAll();
  m_unusedStereoHits.removeAll();
  m_removedHits.removeAll();
  m_circles.removeAll();
  m_tracks.removeAll();
  m_2dTracks.removeAll();
//Liuqg 060917
  for(int i=0;i<19;++i)
    m_unusedAxialHitsOnEachLayer[i].removeAll();
  for(int i=0;i<24;++i)
    m_unusedStereoHitsOnEachLayer[i].removeAll();
  for(int i=0;i<5;++i)
    m_unusedAxialHitsOnEachSuperLayer[i].removeAll();
  for(int i=0;i<6;++i)
    m_unusedStereoHitsOnEachSuperLayer[i].removeAll();
  m_hitsOnInnerSuperLayer.removeAll();
}
 
//
// ... doit .. Section ... This is a main section.
//
int
TCurlFinder::doit(const AList<TMDCWireHit> & axialHits,
          const AList<TMDCWireHit> & stereoHits,
          AList<TTrack> & tracks,
          AList<TTrack> & tracks2D) {
#if (DEBUG_CURL_DUMP+DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
  Belle_event_Manager &evtMgr = Belle_event_Manager::get_manager();
  debugMcFlag = 1;
  if(evtMgr.count() != 0 &&
     evtMgr[0].ExpMC() != 2)debugMcFlag = 0;// not MC
  m_debugCdcFrame = false;
#endif
#if !(DEBUG_CURL_MC)
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)Plot Menu : All Off = 0, Interactive = 1, All On = 2" << std::endl;
  cin >> m_debugPlotFlag;
#endif
  // sub main functions #1, #2, #3
  //...#1
  makeWireHitsListsSegments(axialHits, stereoHits);
#if DEBUG_CURL_SEGMENT
  std::cout << "(TCurlFinder)# of segment = " << m_segmentList.length() << std::endl;
  debugCheckSegments0();
  debugCheckSegments1();
  debugCheckSegments2();
#endif
  //...#2
  if(checkSortSegments() == 0)return 0;
#if DEBUG_CURL_DUMP
  if(m_debugPlotFlag){
    int noPlot = 1;
    if(m_debugPlotFlag == 1){
      std::cout << "(TCurlFinder) Do you want to see Segment Plot? : yes = 1, no = other #" << std::endl;
      cin >> noPlot;
    }
    if(noPlot == 1){
      for(int i=0;i<m_segmentList.length();++i)
    plotSegment(m_segmentList[i]->list(),0);
    }
  }
#endif
  //...#3
  makeCurlTracks(tracks,tracks2D);
#else
  makeWithMC(axialHits, stereoHits, tracks);
#endif
 
  //...iw 2001/01/26...
  unsigned n = tracks2D.length();
  for (unsigned i = 0; i < n; i++)
      tracks2D[i]->quality(TrackQuality2D);
  //...iw end...
  return 0;
}
 
//
// ... Sub Section #1 -- makes segments --
// ... fuction name = "TMakeUnusedWireHitsList"
//
void
TCurlFinder::makeWireHitsListsSegments(const AList<TMDCWireHit> &axialList, 
                       const AList<TMDCWireHit> &stereoList) {  
  // A sub main function ... called by "doit".
  // #0 makes lists.
  // #1 makes segments.
  
  //...makes original lists(axial and stereo)
  //
  //......axial
  unsigned size = axialList.length();
  for(unsigned i=0;i<size;++i){
    if(axialList[i]->reccdc()->tdc > 500) continue; //Liuqg, tmp
    if(axialList[i]->state() & WireHitFindingValid){
      m_allAxialHitsOriginal.append(new TMLink(0,axialList[i]));
      if(axialList[i]->wire()->superLayerId() <= 3)  //origin is 2, Liuqg
    m_hitsOnInnerSuperLayer.append(axialList[i]);
 
    }
  }
  size = m_allAxialHitsOriginal.length();
  for(unsigned i=0;i<size;++i){
    if(!(m_allAxialHitsOriginal[i]->hit()->state() & WireHitUsed)){
      if(m_allAxialHitsOriginal[i]->hit()->state() & WireHitInvalidForFit){
    unsigned newState = m_allAxialHitsOriginal[i]->hit()->state()&(~WireHitInvalidForFit);
    m_allAxialHitsOriginal[i]->hit()->state(newState);
      }
      m_unusedAxialHitsOriginal.append(m_allAxialHitsOriginal[i]);
    }
  }
  m_unusedAxialHits = m_unusedAxialHitsOriginal;
  //......stereo
  size = stereoList.length();
  for(unsigned i=0;i<size;++i){
    if (stereoList[i]->reccdc()->tdc > 500) continue; //Liuqg, tmp to exclude looping tracks
    if(stereoList[i]->state() & WireHitFindingValid){
      m_allStereoHitsOriginal.append(new TMLink(0,stereoList[i]));
      if(stereoList[i]->wire()->superLayerId() <= 3)  //origin is 2, Liuqg
    m_hitsOnInnerSuperLayer.append(stereoList[i]);
    }
  }
  size = m_allStereoHitsOriginal.length();
  for(unsigned i=0;i<size;++i){
    if(!(m_allStereoHitsOriginal[i]->hit()->state() & WireHitUsed)){
      if(m_allStereoHitsOriginal[i]->hit()->state() & WireHitInvalidForFit){
    unsigned newState = m_allStereoHitsOriginal[i]->hit()->state()&(~WireHitInvalidForFit);
    m_allStereoHitsOriginal[i]->hit()->state(newState);
      }
      m_unusedStereoHitsOriginal.append(m_allStereoHitsOriginal[i]);
    }
  }
  m_unusedStereoHits = m_unusedStereoHitsOriginal;
 
  //...shares unsed hit wires to each layer.
  size = m_unusedAxialHitsOriginal.length();
  for(unsigned i=0;i<size;++i){
    m_unusedAxialHitsOnEachLayer[m_unusedAxialHitsOriginal[i]->hit()->wire()->
                axialStereoLayerId()].append(m_unusedAxialHitsOriginal[i]);
  }
  size = m_unusedStereoHitsOriginal.length();
  for(unsigned i=0;i<size;++i){
    m_unusedStereoHitsOnEachLayer[m_unusedStereoHitsOriginal[i]->hit()->wire()->
                 axialStereoLayerId()].append(m_unusedStereoHitsOriginal[i]);
  }
  
  //...sets pointers to neighboring hit wires of each TMLink
  linkNeighboringWires(m_unusedAxialHitsOnEachLayer,19);   //Belle 32
  linkNeighboringWires(m_unusedStereoHitsOnEachLayer,24);  //Belle 18
  
  //...makes _unusedSuperAxialHits and _unusedSuperStereoHits
  createSuperLayer();
  //...makes segments by linking neighboring hit wires in the same super layer
  m_segmentList.removeAll();
  for(unsigned i=0;i<5;++i) //Belle 6
    if(m_unusedAxialHitsOnEachSuperLayer[i].length() > 0)
      createSegments(m_unusedAxialHitsOnEachSuperLayer[i]);
  for(unsigned i=0;i<6;++i) //Belle 5
    if(m_unusedStereoHitsOnEachSuperLayer[i].length() > 0)
      createSegments(m_unusedStereoHitsOnEachSuperLayer[i]);
 
}
 
void
TCurlFinder::linkNeighboringWires(AList<TMLink> *list, const unsigned num) {
  // Axial(num == 19) and Stereo(num == 24).
  // ...sets pointers to neighboring wires 
  // in "neighbor" of "AList<TMLink> *list" element
  
  for(int i=0;i<num;++i){
    if(list[i].length() == 0)continue;
    for(int j=0;j<list[i].length();++j){
      //find two links of list[i][j]
      //...inner layer
      if(num == 19){
        if( i == 0  || i == 4  || i == 8 ||
        i == 12 || i == 16)goto outer;
      }else if(num == 24){
        if( i == 0 || i == 4 ||
        i == 8 || i == 12 || i == 16 || i == 20 )goto outer;
      }
      
      for(int k=0;k<list[i-1].length();++k){
        if(list[i-1][k]->hit()->wire()->localId() ==
           list[i][j]->wire()->neighbor(1)->localId())
          setNeighboringWires(list[i][j], list[i-1][k]);
        else if(list[i-1][k]->hit()->wire()->localId() ==
             list[i][j]->wire()->neighbor(0)->localId())
            setNeighboringWires(list[i][j], list[i-1][k]);
      }//k
    outer:
      //...outer layer
      if(num == 19){
        if( i ==  3 || i == 7 || i == 11 ||
        i == 15 || i == 18)goto same;
      }else if(num == 24){
    if( i ==  3 || i == 7 || i == 11 ||
            i == 15 || i == 19 || i == 23 )goto same;
      }
      for(int k=0;k<list[i+1].length();++k){
        if(list[i+1][k]->hit()->wire()->localId() ==
           list[i][j]->hit()->wire()->neighbor(4)->localId())
          setNeighboringWires(list[i][j], list[i+1][k]);
        else if(list[i+1][k]->wire()->localId() ==
             list[i][j]->wire()->neighbor(5)->localId())
            setNeighboringWires(list[i][j], list[i+1][k]);
      }//k
    same:
      //...same layer
      for(int k=0;k<list[i].length();++k){
        if(list[i][k]->hit()->wire()->localId() ==
           list[i][j]->hit()->wire()->localIdForPlus()+1){
      setNeighboringWires(list[i][j], list[i][k]);
        }else if(list[i][k]->hit()->wire()->localId() ==
         list[i][j]->hit()->wire()->localIdForMinus()-1){
      setNeighboringWires(list[i][j], list[i][k]);
        }
      }//k
    }//j
  }//i
}
 
void
TCurlFinder::setNeighboringWires(TMLink *list, const TMLink *next) {
  // ...sets a neighboring wire of "list".
  // Its candidate is "next".
  for(int i=0;i<6;++i){
    if(!list->neighbor(i)){
      list->neighbor(i,const_cast<TMLink*>(next));
      break;
    }
  }
}
 
void
TCurlFinder::createSuperLayer(void) {
//Liuqg
  for(int i=0;i<4;++i){
    if(m_unusedAxialHitsOnEachLayer[i].length() > 0)
      m_unusedAxialHitsOnEachSuperLayer[0].append(m_unusedAxialHitsOnEachLayer[i]);
    if(m_unusedAxialHitsOnEachLayer[i+4].length() > 0)
      m_unusedAxialHitsOnEachSuperLayer[1].append(m_unusedAxialHitsOnEachLayer[i+4]);
    if(m_unusedAxialHitsOnEachLayer[i+8].length() > 0)
      m_unusedAxialHitsOnEachSuperLayer[2].append(m_unusedAxialHitsOnEachLayer[i+8]);
    if(m_unusedAxialHitsOnEachLayer[i+12].length() > 0)
      m_unusedAxialHitsOnEachSuperLayer[3].append(m_unusedAxialHitsOnEachLayer[i+12]);
    if(m_unusedAxialHitsOnEachLayer[i+16].length() > 0 && i+16 < 19)
      m_unusedAxialHitsOnEachSuperLayer[4].append(m_unusedAxialHitsOnEachLayer[i+16]);
  }
  for(int i=0;i<4;++i){
    if(m_unusedStereoHitsOnEachLayer[i].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[0].append(m_unusedStereoHitsOnEachLayer[i]);
    if(m_unusedStereoHitsOnEachLayer[i+4].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[1].append(m_unusedStereoHitsOnEachLayer[i+4]);
    if(m_unusedStereoHitsOnEachLayer[i+8].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[2].append(m_unusedStereoHitsOnEachLayer[i+8]);
    if(m_unusedStereoHitsOnEachLayer[i+12].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[3].append(m_unusedStereoHitsOnEachLayer[i+12]);
    if(m_unusedStereoHitsOnEachLayer[i+16].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[4].append(m_unusedStereoHitsOnEachLayer[i+16]);
    if(m_unusedStereoHitsOnEachLayer[i+20].length() > 0)
      m_unusedStereoHitsOnEachSuperLayer[5].append(m_unusedStereoHitsOnEachLayer[i+20]);
  }
}
 
void
TCurlFinder::createSegments(AList<TMLink> &list) {
  // ...makes segments from AList<TMLink> &list, in every superlayer. 
  // These segments are add to _segmentList.(# of segments >= MIN_SEGMENT)
  AList<TMLink> seedStock;
  do{
    TSegmentCurl *segment = new TSegmentCurl(list[0]->hit()->wire()->superLayerId(),
                         maxLocalLayerId(list[0]->hit()->wire()
                                 ->superLayerId()));
 
    segment->append(list[0]);    
    TMLink *seed = list[0];
    list.remove(seed);
  next:
    searchSegment(seed, list, seedStock, segment);
    if(seedStock.length() > 0){
      seed = seedStock[0];
      seedStock.remove(seed);
      goto next;
    }else if(segment->size() >= m_param.MIN_SEGMENT){
      segment->update();
      m_segmentList.append(segment);
#if DEBUG_CURL_DUMP
      std::cout << "Segment # = " << m_segmentList.length() << std::endl;
      segment->dump();
#endif
    }else{
      delete segment;
    }
  }while(list.length() > 0);
}
 
void
TCurlFinder::searchSegment(TMLink *seed, AList<TMLink> &list, 
               AList<TMLink> &seedStock, TSegmentCurl *segment) {
  for(int i=0;i<6;++i){
    if(seed->neighbor(i)){
      if(!findLink(seed->neighbor(i),list))continue;
      segment->append(seed->neighbor(i));
      seedStock.append(seed->neighbor(i));
      list.remove(seed->neighbor(i));
    }else{
      break;
    }
  }
}
 
TMLink *
TCurlFinder::findLink(const TMLink *seed, const AList<TMLink> &list) {
  // This is to search "TMLink *seed" in "AList<TMLink> &list".
  // Return is when found, the "TMLink *list[i]"
  // when not found, NULL.
  unsigned size = list.length();
  if(size == 0)return NULL;
  for(unsigned i=0;i<size;++i){
    if(seed == list[i])return list[i];
  }
  return NULL;
}
 
//
// ... Sub Section #2 -- checks and sorts segments --
// ... fuction name = "checkSortSegments"
//
int
TCurlFinder::checkSortSegments(void) {
  // A sub main function ...called by "doit".
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)checking and sorting segments..." << std::endl;
#endif
  unsigned length = m_segmentList.length();
  if(length == 0)return 0;
  checkExceptionalSegmentsType03();//...exception #3
  //checkExceptionalSegmentsType02();//...exception #2
  checkExceptionalSegmentsType01();//...exception #1
  m_segmentList.sort(sortBySequentialLength);
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)...done check and sort of segments." << std::endl;
#endif
  return 1;
}
 
void
TCurlFinder::checkExceptionalSegmentsType03(void) {
  int max = m_param.MAX_FULLWIRE;
  int nMinWires;
  if(max == 7)nMinWires = 21;
  else if(max == 6)nMinWires = 19;
  else if(max == 5)nMinWires = 18;
  else if(max == 4)nMinWires = 16;
  else if(max == 3)nMinWires = 14;
  else if(max == 2)nMinWires = 12;
  else if(max == 1)nMinWires = 10;
  else if(max == 0)nMinWires = 7;
 
  AList<TSegmentCurl> removeList;
  for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
    if(m_segmentList[i]->size() >= nMinWires){
      unsigned nWires = m_segmentList[i]->size();
      unsigned n6Wires = 0;
      for(unsigned j=0;j<nWires;++j){
    if(((m_segmentList[i]->list())[j])->neighbor(5))++n6Wires;
    if(n6Wires > max)break;
      }
      if(n6Wires <= max)continue;
      removeList.append(m_segmentList[i]);
#if DEBUG_CURL_SEGMENT
      writeSegment(m_segmentList[i]->list(),3);
#endif
    }
  }
  if(removeList.length() >= 1){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)removing large segments: # = " << removeList.length() << std::endl;
#endif
    m_segmentList.remove(removeList);
    HepAListDeleteAll(removeList);    
  }
}
 
void
TCurlFinder::checkExceptionalSegmentsType02(void) {
  int max  = 10;
  int hmax = 5;
  AList<TSegmentCurl> removeList;
  for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
    int lSize = max*3+hmax*3;
    int lNum  = 3;
    if(m_segmentList[i]->superLayerId() == 1 ||
       m_segmentList[i]->superLayerId() == 3){
      lSize = max*2+hmax;
      lNum = 2;
    }
    if(m_segmentList[i]->superLayerId() == 5 ||
       m_segmentList[i]->superLayerId() == 7 ||
       m_segmentList[i]->superLayerId() == 9){
      lSize = max*2+hmax*2;
      lNum  = 2;
    }
    if(m_segmentList[i]->superLayerId() == 4 ||
       m_segmentList[i]->superLayerId() == 6 ||
       m_segmentList[i]->superLayerId() == 8 ||
       m_segmentList[i]->superLayerId() == 10)lSize = max*3+hmax*2;
    if(m_segmentList[i]->size() < lSize)continue;
    int nL = 0;
    for(unsigned j=0,size=m_segmentList[i]->maxLocalLayerId();j<size;++j){
      if(m_segmentList[i]->sizeOfLayer(j) >= max)++nL;
    }
    if(nL < lNum)continue;
    removeList.append(m_segmentList[i]);
    //plotSegment(m_segmentList[i]->list(),0);
#if DEBUG_CURL_SEGMENT
    //writeSegment(m_segmentList[i]->list(),2);
#endif
  }
  if(removeList.length() >= 1){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)removing large segments: # = " << removeList.length() << std::endl;
#endif
    m_segmentList.remove(removeList);
    HepAListDeleteAll(removeList);    
  }
}
 
void
TCurlFinder::checkExceptionalSegmentsType01(void) {
  for(unsigned i = 0, length = m_segmentList.length(); i < length; ++i){
    if(m_segmentList[i]->maxLocalLayerId() != m_segmentList[i]->layerIdOfMaxSeq() &&
       m_segmentList[i]->maxSeq() >= m_param.MIN_SEQUENCE){
      unsigned innerHits = 0;
      if(m_segmentList[i]->layerIdOfMaxSeq() == 0)continue;
      TSegmentCurl *outer = new TSegmentCurl(m_segmentList[i]->superLayerId(), 
                         m_segmentList[i]->maxLocalLayerId());
      for(unsigned j = 0, size = m_segmentList[i]->size(); j < size; ++j){
    if(m_segmentList[i]->layerIdOfMaxSeq()+1 <= 
       (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId() &&
       (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId() <= 
       m_segmentList[i]->maxLocalLayerId()){
      outer->append((m_segmentList[i]->list())[j]);
    }else if(m_segmentList[i]->layerIdOfMaxSeq()-1 >= 
         (m_segmentList[i]->list())[j]->hit()->wire()->localLayerId()){
      ++innerHits;
    }
      }
      if(innerHits != 0 && outer->size() != 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)removing some wires in the segment." << std::endl;
#endif
#if DEBUG_CURL_SEGMENT
    //writeSegment(m_segmentList[i]->list(),1);
#endif
        m_segmentList[i]->remove(const_cast< AList<TMLink>& >(outer->list()));
    outer->removeAll();
    delete outer;
      }else{
    outer->removeAll();
    delete outer;
      }
    }
  }
}
 
//
// ... Sub Section #3 -- makes curl tracks --
// ... fuction name = "makeCurlTracks"
//
 
//............................................
//..............2D+3D Manage Parts............
//............................................
 
void
TCurlFinder::makeCurlTracks(AList<TTrack> &tracks,
                AList<TTrack> &tracks2D) {
  AList<TSegmentCurl> segmentList = m_segmentList;
  //cout<<"segments: "<<m_segmentList.length()<<endl;
//  if(m_param.SVD_RECONSTRUCTION)m_builder.setSvdClusters();
  for(unsigned i = 0, size = m_segmentList.length(); i < size; ++i){
    TCircle *circle = make2DTrack(segmentList[i]->list(), segmentList, 1);   //order by sequential num
    if(circle){
      AList<TMLink> tmp(circle->links());
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)  2D:Created Circle!!!" << std::endl;
      if(m_debugPlotFlag){
    int noPlot = 1;
    if(m_debugPlotFlag == 1){
      std::cout << "(TCurlFinder)   Do you want to see Circle Plot(2D)? : yes = 1, no = other #" << std::endl;
      cin >> noPlot;
    }
    if(noPlot == 1)plotCircle(*circle,0);
      }
#endif
 
/*      if (tmp.length()>0) {
    cout<<"circle ok.............."<<endl;
    for(int j = 0; j < tmp.length(); ++j){
      TMLink *ll = tmp[j];
      cout<<"layerId: "<<ll->wire()->layerId()
        <<"  localId: "<<ll->wire()->localId()<<endl;
    }
      }
*/
      if(TCircle *dividedCircle = dividing2DTrack(circle)){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)   2D:dividing...good...2 Circles!!!" << std::endl;
#endif
    TTrack *track1(NULL), *track2(NULL);
    int ok2d[2] = { 0, 0 };
    int ok3d[2] = { 0, 0 };
    //cout<<"trac2D: "<<trace2DTrack(circle)<<"  check2D: "<< check2DCircle(circle)
      //<<"  fit: "<< circle->fitForCurl(1)<<endl;
    if(trace2DTrack(circle) &&
       check2DCircle(circle) &&
       circle->fitForCurl(1) != -1){
      ok2d[0] = 1;
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)    2D:Success Circle Fit!!!" << std::endl;
#endif
      track1 = make3DTrack(circle, segmentList);
    }
    //cout<<"trac2D: "<<trace2DTrack(dividedCircle)<<"  check2D: "<< check2DCircle(dividedCircle)
      //<<"  fit: "<< dividedCircle->fitForCurl(1)<<endl;
    if(trace2DTrack(dividedCircle) &&
       check2DCircle(dividedCircle) &&
       dividedCircle->fitForCurl(1) != -1){
      ok2d[1] = 1;
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)    2D:Success Circle Fit!!!" << std::endl;
#endif
      track2 = make3DTrack(dividedCircle, segmentList);
    }
    if(track1 && track2){
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)      3D:Create Track!!! in track1 && track2" << std::endl;
#endif
      salvage3DTrack(track1,true);
      salvage3DTrack(track2,true);
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)      1:dz = " << track1->helix().dz() 
           << ", 2:dz = " << track2->helix().dz() << std::endl;
#endif
      if(m_param.DETERMINE_ONE_TRACK){
        if(fabs(track1->helix().dz()) < fabs(track2->helix().dz())){
          m_tracks.remove(track2);
          if(merge3DTrack(track1, tracks))
        if(check3DTrack(track1) &&
           trace3DTrack(track1)){
          ok3d[0] = 1;
          ok3d[1] = 1;
          mask3DTrack(track1,tmp);
          tmp.append(track1->links());
        }
        }else{        
          m_tracks.remove(track1);
          if(merge3DTrack(track2, tracks))
        if(check3DTrack(track2) &&
           trace3DTrack(track2)){
          ok3d[0] = 1;
          ok3d[1] = 1;
          mask3DTrack(track2,tmp);
          tmp.append(track2->links());
        }
        }
      }else{
        int isSaved[2] = { 0, 0 };
        if(merge3DTrack(track1, tracks)){
          if(check3DTrack(track1) &&
         trace3DTrack(track1)){
        ok3d[0] = 1;
        mask3DTrack(track1,tmp);
        tmp.append(track1->links());
        isSaved[0] = 1;
          }
        }
        if(merge3DTrack(track2, tracks)){
          if(check3DTrack(track2) &&
         trace3DTrack(track2)){
        ok3d[1] = 1;
        mask3DTrack(track2,tmp);
        tmp.append(track2->links());
        isSaved[1] = 1;
          }
        }
        if(isSaved[0] == 1 && isSaved[1] == 1){
          track1->daughter(track2);
          track2->daughter(track1);
        }
      }
    }else if(track1){
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)      3D:Create Track!!! in track1" << std::endl;
#endif
      salvage3DTrack(track1,true);
      if(merge3DTrack(track1, tracks))
        if(check3DTrack(track1) &&
           trace3DTrack(track1)){
          ok3d[0] = 1;
          mask3DTrack(track1,tmp);
          tmp.append(track1->links());
        }
    }else if(track2){
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)      3D:Create Track!!! in track2" << std::endl;
#endif
      salvage3DTrack(track2,true);
      if(merge3DTrack(track2, tracks))
        if(check3DTrack(track2) &&
           trace3DTrack(track2)){
          ok3d[1] = 1;
          mask3DTrack(track2,tmp);
          tmp.append(track2->links());
        }
    }
 
    if(m_param.OUTPUT_2DTRACKS){
      // When 2d is OK but 3d is BAD, a 2dtrk is saved.
      if(ok2d[0] == 1 && ok3d[0] == 0){
        removeStereo(*circle);
        double chi2_2d;
        int    ndf_2d;
        if(fitWDD(*circle,chi2_2d,ndf_2d)){
          TTrack *trk2d = new TTrack(*circle);
          trk2d->_ndf   = ndf_2d;
          trk2d->_chi2  = chi2_2d;
          m_2dTracks.append(trk2d);
          m_allTracks.append(trk2d);
        }
#if DEBUG_CURL_DUMP
        else{
          std::cout << "(TCurlFinder)   2D:fit with drift information!!!" << std::endl;
        }
#endif
      }
      if(ok2d[1] == 1 && ok3d[1] == 0){
        removeStereo(*dividedCircle);
        double chi2_2d;
        int    ndf_2d;
        if(fitWDD(*dividedCircle,chi2_2d,ndf_2d)){
          TTrack *trk2d = new TTrack(*dividedCircle);
          trk2d->_ndf   = ndf_2d;
          trk2d->_chi2  = chi2_2d;
          m_2dTracks.append(trk2d);
          m_allTracks.append(trk2d);
        }
#if DEBUG_CURL_DUMP
        else{
          std::cout << "(TCurlFinder)   2D:fit with drift information!!!" << std::endl;
        }
#endif
      }
    }
 
      }else{
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)   2D:dividing...no good...1 Circles!!!" << std::endl;
#endif
    int ok2d = 0;
    int ok3d = 0;
    //cout<<"trac2D: "<<trace2DTrack(circle)<<"  check2D: "<< check2DCircle(circle)
      //<<"  fit: "<< circle->fitForCurl(1)<<endl;
    if(trace2DTrack(circle) &&
       check2DCircle(circle) &&
       circle->fitForCurl(1) != -1){
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)    2D:Success Circle Fit!!!" << std::endl;
#endif
      ok2d = 1;
      TTrack *track3 = make3DTrack(circle, segmentList);
      if(track3){
#if DEBUG_CURL_DUMP
        std::cout << "(TCurlFinder)      3D:Create Track!!! in track3" << std::endl;
#endif
        salvage3DTrack(track3,true);
        if(merge3DTrack(track3, tracks))
          if(check3DTrack(track3) &&
         trace3DTrack(track3)){
        ok3d = 1;
        mask3DTrack(track3,tmp);
        tmp.append(track3->links());
          }
      }
 
      if(m_param.OUTPUT_2DTRACKS){    
        // When 2d is OK but 3d is BAD, a 2dtrk is saved.
        if(ok2d == 1 && ok3d == 0){
          removeStereo(*circle);
          double chi2_2d;
          int    ndf_2d;
          if(fitWDD(*circle,chi2_2d,ndf_2d)){
        TTrack *trk2d = new TTrack(*circle);
        trk2d->_ndf   = ndf_2d;
        trk2d->_chi2  = chi2_2d;
        m_2dTracks.append(trk2d);
        m_allTracks.append(trk2d);
          }
#if DEBUG_CURL_DUMP
          else{
        std::cout << "(TCurlFinder)   2D:fit with drift information!!!" << std::endl;
          }
#endif
        }
      }
    }
      }
      m_unusedAxialHits.remove(tmp);
      m_unusedStereoHits.remove(tmp);
      for(unsigned ii=0, nsize=m_segmentList.length();ii<nsize;++ii){
    m_segmentList[ii]->remove(tmp);
    if(m_segmentList[ii]->list().length() < m_param.MIN_SEGMENT)
      m_segmentList[ii]->removeAll();
      }
    }
    segmentList[i]->removeAll();
  }
 
  // Check 2D trk's wires
  //check2DTracks(); // ... not need because of the current algorithm
 
  // 3D & 2D trks information
  assignTracks();
 
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)MDC Rec Track # 3D = " << m_tracks.length()
        << ", 2D = " << m_2dTracks.length() << std::endl;
  std::cout << "3D Track List" << std::endl;
  for(int j=0;j<m_tracks.length();++j){
      m_tracks[j]->setFinderType(3);
    unsigned nA = 0, nS = 0;
    unsigned nAOK = 0, nSOK = 0;
    for(unsigned i=0,size=m_tracks[j]->nLinks();i<size;++i){
      if(m_tracks[j]->links()[i]->wire()->stereo())++nS;
      else ++nA;
      if(/*!(m_tracks[j]->links()[i]->hit()->state() & WireHitInvalidForFit) &&*/
     (m_tracks[j]->links()[i]->hit()->state() & WireHitFittingValid)){
    if(m_tracks[j]->links()[i]->wire()->stereo())++nSOK;
    else ++nAOK;
      }
    }
    std::cout << "(TCurlFinder)  #" << j << ": wire info...A+S: " << m_tracks[j]->nLinks() 
          << ", A: " << nAOK << "/" << nA
          << ", S: " << nSOK << "/" << nS << std::endl;
    if(m_tracks[j]->daughter())
      std::cout << "(TCurlFinder)  Relation = EXIST" << std::endl;
    else
      std::cout << "(TCurlFinder)  Relation = NO EXIST" << std::endl;
  }
  std::cout << "2D Track List" << std::endl;
  for(int j=0;j<m_2dTracks.length();++j){
    unsigned nA = 0, nS = 0;
    unsigned nAOK = 0, nSOK = 0;
    for(unsigned i=0,size=m_2dTracks[j]->nLinks();i<size;++i){
      if(m_2dTracks[j]->links()[i]->wire()->stereo())++nS;
      else ++nA;
      if(/*!(m_2dTracks[j]->links()[i]->hit()->state() & WireHitInvalidForFit) &&*/
         (m_2dTracks[j]->links()[i]->hit()->state() & WireHitFittingValid)){
        if(m_2dTracks[j]->links()[i]->wire()->stereo())++nSOK;
        else ++nAOK;
      }
    }
    std::cout << "(TCurlFinder)  #" << j << ": wire info...A+S: " << m_2dTracks[j]->nLinks() 
          << ", A: " << nAOK << "/" << nA
          << ", S: " << nSOK << "/" << nS 
          << ", Chi2: " << m_2dTracks[j]->chi2()
          << ", Ndf: " << m_2dTracks[j]->ndf() << std::endl;
    if(m_2dTracks[j]->daughter())
      std::cout << "(TCurlFinder)  Relation = EXIST" << std::endl;
    else
      std::cout << "(TCurlFinder)  Relation = NO EXIST" << std::endl;
  }
#endif
 
  // 3D trks
  m_allTracks.remove(m_tracks);
  checkRelation(m_tracks);
  tracks.append(m_tracks);
  if(m_param.OUTPUT_2DTRACKS){
    // 2D trks
    m_allTracks.remove(m_2dTracks);
    tracks2D.append(m_2dTracks);
  }
}
 
void
TCurlFinder::check2DTracks(void)
{
  if(m_2dTracks.length() == 0)return;
  AList<TMLink> allWires_3Dtrks;
  for(int i=0;i<m_tracks.length();++i){
    allWires_3Dtrks.append(m_tracks[i]->links());
  }
  //cout << "ALL Wire(3D) " << allWires_3Dtrks.length() << endl;
 
  for(int i=0;i<m_2dTracks.length();++i){
    AList<TMLink> usedWires;
    for(int j=0;j<m_2dTracks[i]->nLinks();++j){
      int ok = 1;
      for(int k=0;k<allWires_3Dtrks.length();++k){
    if(m_2dTracks[i]->links()[j]->wire()->id() ==
       allWires_3Dtrks[k]->wire()->id()){
      ok = 0;
      break;
    }
      }
      if(ok == 0){
    usedWires.append(m_2dTracks[i]->links()[j]);
      }
    }
    //cout << i << " : used # " << usedWires.length() 
    // << ", all # " << m_2dTracks[i]->nLinks() << endl;
    m_2dTracks[i]->remove(usedWires);
  }
}
 
void
TCurlFinder::checkRelation(AList<TTrack> &list) {
  unsigned nT = list.length();
  if(nT <= 1)return;
  for(unsigned i=0;i<nT;++i){
    if(list[i]->daughter()){
      int isHere = 0;
      for(unsigned j=0;j<nT;++j){
    if(i != j &&
       list[i]->daughter() == list[j]){
      isHere = 1;
      break;
    }
      }
      if(isHere == 0){
    list[i]->daughter(NULL);
      }
    }
  }
}
 
TCircle *
TCurlFinder::dividing2DTrack(TCircle *circle) {
  AList<TMLink> positive, negative;
  for(unsigned i = 0, size = circle->nLinks(); i < size; ++i){
    if(circle->center().x()*circle->links()[i]->hit()->wire()->xyPosition().y() -
       circle->center().y()*circle->links()[i]->hit()->wire()->xyPosition().x() > 0.){
      positive.append(circle->links()[i]);
    }else{
      negative.append(circle->links()[i]);
    }
  }
  if(positive.length() > negative.length()){
    circle->remove(negative);
    circle->property(1.,fabs(circle->radius()),circle->center());
    if(negative.length() >= 3){
      TCircle *new_circle = new TCircle(negative);
      m_allCircles.append(new_circle);
      new_circle->property(-1.,-1.*fabs(circle->radius()),circle->center());
      return new_circle;
    }else{
      return NULL;
    }
  }else{
    circle->remove(positive);
    circle->property(-1.,-1.*fabs(circle->radius()),circle->center());
    if(positive.length() >= 3){
      TCircle *new_circle = new TCircle(positive);
      m_allCircles.append(new_circle);
      new_circle->property(1.,fabs(circle->radius()),circle->center());
      return new_circle;
    }else{
      return NULL;
    }
  }
}
 
void
TCurlFinder::assignTracks(void) {
  // 3D trks
  for(int i=0,size=m_tracks.length();i<size;++i) {
    m_tracks[i]->assign(WireHitCurlFinder);
    m_tracks[i]->finder(TrackCurlFinder);
//    m_tracks[i]->assign(WireHitCurlFinder, TrackCurlFinder | TrackValid | Track3D);
  }
 
  // 2D trks
  for(int i=0,size=m_2dTracks.length();i<size;++i) {
    m_2dTracks[i]->assign(WireHitCurlFinder);
    m_2dTracks[i]->finder(TrackCurlFinder);
  }
}
 
extern "C" int 
TCurlFinder_doubleCompare(const void *i, const void *j){
  if(*(static_cast<const double*>(i)) > *(static_cast<const double*>(j)))return 1;
  if(*(static_cast<const double*>(i)) < *(static_cast<const double*>(j)))return -1;
  return 0;
}
 
int
TCurlFinder::trace2DTrack(TCircle *circle) {
  // 0 is bad, 1 is good
  unsigned nSize = circle->links().length();
  if(nSize == 0)return 0;
  double r = fabs(circle->radius());
  if(r < 0.01)return 0; // to reject "r=0cm" circles.
  double cx = circle->center().x();
  double cy = circle->center().y();
  double th = atan2(-cy,-cx);
  if(th < 0.)th += 2.*M_PI;
 
  unsigned innerOK = 0;
  double *angle = new double [circle->links().length()];
  for(unsigned i=0,size=nSize;i<size;++i){
    double th_r = atan2(circle->links()[i]->wire()->xyPosition().y()-cy,
            circle->links()[i]->wire()->xyPosition().x()-cx);
    if(th_r < 0.)th_r += 2.*M_PI;
    double diff = th_r-th+2.*M_PI;
    if(th_r > th)diff = th_r-th;
    if(circle->links()[i]->wire()->superLayerId() <= 2)innerOK = 1;
    angle[i] = diff;
  }
  qsort(angle,nSize,sizeof(double),TCurlFinder_doubleCompare);//chiisai-jun
  double maxDiffAngle = 0.;
  unsigned maxIndex = 0;
  for(unsigned i=0,size=nSize;i<size-1;++i){
    if(angle[i+1]-angle[i] > maxDiffAngle){
      maxDiffAngle = angle[i+1]-angle[i];
      maxIndex = i;
    }
  }
  delete [] angle;
  //std::cout << "2D TRACE : maxDifAngle = " << maxDiffAngle 
    // << ", r = " << r << ", dist = " << r*maxDiffAngle << std::endl;
 
//  if(r*maxDiffAngle > m_param.TRACE2D_DISTANCE)return 0;  //Liuqg, tmp. BES layer distribution is not uniform.
  if(innerOK == 1)return 1;
 
  double q = circle->radius() > 0. ? 1. : -1;
  for(unsigned i=0,size=m_hitsOnInnerSuperLayer.length();i<size;++i){
    double mag = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx,
              m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy);
    if(fabs(mag-r) < m_param.TRACE2D_FIRST_SUPERLAYER &&
       q*(cx*m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-
      cy*m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()) > 0.){
      return 1;
    }
  }
  return 0;
}
 
bool
TCurlFinder::check2DCircle(TCircle *circle) {
  unsigned nA(nAxialHits(fabs(circle->radius())*2.0));
  
  unsigned nMA = static_cast<unsigned>(floor(m_param.RATIO_USED_WIRE*static_cast<double>(nA)));
  if(nMA < 3)nMA = 3;
 
  unsigned nAhits(0), nShits(0);
  for(unsigned i=0,size=circle->nLinks();i<size;++i){
    if((circle->links()[i])->wire()->axial())++nAhits;
    else ++nShits;
  }
#if DEBUG_CURL_DUMP
  if(nAhits < nMA){
    std::cout << "(TCurlFinder)    2D:Fail...checking axial wires # = " 
     << nAhits << " < " << nMA << std::endl;
  }
#endif
  if(nAhits >= nMA)return true;
  return false;
}
 
bool
TCurlFinder::check3DTrack(TTrack *track) {
  trace3DTrack(track);
  unsigned nA = 0, nS = 0;
  for(unsigned i=0,size=track->nLinks();i<size;++i){
    if(!(track->links()[i]->hit()->state() & WireHitFittingValid))continue;
    if(track->links()[i]->wire()->stereo())++nS;
    else ++nA;
    if(nA >= 3 && nS >= 2)return true;
  }
  m_tracks.remove(track);
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)       3D:Checked...Fail...removing this track. Valid Axial # = " 
       << nA << ", Stereo # = " << nS << std::endl;
#endif
  return false;
}
 
int
TCurlFinder::trace3DTrack(TTrack *track) {
  // 0 is bad, 1 is good
  unsigned nSize = track->links().length();
  if(nSize == 0){
    m_tracks.remove(track);
    return 0;
  }
  double r = fabs(track->helix().radius());
  if(r < 0.01){
    m_tracks.remove(track);
    return 0; // to reject "r=0cm" circles.
  }
  double cx = track->helix().center().x();
  double cy = track->helix().center().y();
  double th = atan2(-cy,-cx);
  if(th < 0.)th += 2.*M_PI;
 
  double *angle = new double [track->links().length()];
  for(unsigned i=0,size=nSize;i<size;++i){
    double th_r = atan2(track->links()[i]->positionOnTrack().y()-cy,
            track->links()[i]->positionOnTrack().x()-cx);
    if(th_r < 0.)th_r += 2.*M_PI;
    double diff = th_r-th+2.*M_PI;
    if(th_r > th)diff = th_r-th;
    angle[i] = diff;
  }
  qsort(angle,nSize,sizeof(double),TCurlFinder_doubleCompare);//chiisai-jun
  double maxDiffAngle = 0.;
  unsigned maxIndex = 0;
  for(unsigned i=0,size=nSize;i<size-1;++i){
    if(angle[i+1]-angle[i] > maxDiffAngle){
      maxDiffAngle = angle[i+1]-angle[i];
      maxIndex = i;
    }
  }
  delete [] angle;
  /* std::cout << "3D TRACE : maxDifAngle = " << maxDiffAngle 
     << ", r = " << r << ", dist = " << r*maxDiffAngle << std::endl; */
  if(r*maxDiffAngle > m_param.TRACE3D_DISTANCE){
     m_tracks.remove(track);
    return 0;
  }else{
    return 1;
  }
}
 
void
TCurlFinder::mask3DTrack(TTrack *track,
             AList<TMLink> &maskList) {
  double r(fabs(track->helix().radius()));
  double cx(track->helix().center().x());
  double cy(track->helix().center().y());
 
  AList<TMLink> list(m_unusedAxialHits);
  list.append(m_unusedStereoHits);
  list.remove(track->links());
  list.sort(SortByWireId);
 
  AList<TMLink> removeList;
  for(unsigned i=0,size=list.length();i<size;++i){
    double d = distance(*track, *(list[i]));
    if(d < m_param.MASK_DISTANCE){
      HepPoint3D tmp(d, 0., 0.);
      list[i]->position(tmp);
      removeList.append(list[i]);
    }
  }
 
  int pLayerId1 = static_cast<int>(layerId(2.*r));
  if(pLayerId1 != 50)pLayerId1 -= 1;//hard coding parameter
  int pLayerId2 = pLayerId1+2;      //hard coding parameter
 
  AList<TMLink> preCand, cand;
  while(removeList.length()){
    preCand.removeAll();
    preCand.append(removeList[0]);
    if(removeList.length() >= 2){
      for(unsigned j=1,size=removeList.length();j<size;++j){
    if(removeList[0]->wire()->layerId() == removeList[j]->wire()->layerId()){
      for(unsigned k=0,num=preCand.length();k<num;++k){
        if(preCand[k]->wire()->localIdForPlus()+1 == removeList[j]->wire()->localId()){
          preCand.append(removeList[j]);
          break;
        }
      }
    }
      }
#if 1
      // new
      if(preCand[0]->wire()->layerId() >= pLayerId1 && 
     preCand[0]->wire()->layerId() <= pLayerId2){
    cand.append(preCand);
      }else if(preCand.length() == 2){//hard coding parameter
    cand.append(preCand);
      }else if(preCand.length() == 1){
    cand.append(preCand[0]);
      }
#else
      if(preCand.length() == 1){
    if(preCand[0]->position().x() < MASK_DISTANCE)cand.append(preCand[0]);
      }else{
    if(preCand[0]->wire()->layerId() >= pLayerId1 && 
       preCand[0]->wire()->layerId() <= pLayerId2){
      cand.append(preCand);
    }else if(preCand.length() == 2){//hard coding parameter
      cand.append(preCand);
    }
      }
#endif
    }else{
      cand.append(removeList[0]);
    }
    removeList.remove(removeList[0]);
    removeList.remove(cand);
  }
  maskList.append(cand);
}
 
TTrack*
TCurlFinder::merge3DTrack(TTrack *track, AList<TTrack> &confTracks) {
  if(!m_param.MERGE_EXE)return track;
 
  AList<TTrack> tracks(confTracks);
  tracks.append(m_tracks);
  tracks.remove(track);
  if(tracks.length() == 0)return track;
 
  double r  = track->helix().radius();
  double cx = track->helix().center().x();
  double cy = track->helix().center().y();
 
  unsigned bestIndex = 0;
  double   bestDiff = 1.0e+20;
  double   R, cX, cY;
  for(unsigned i=0,size=tracks.length();i<size;++i){
    R  = fabs(tracks[i]->helix().radius());
    cX = tracks[i]->helix().center().x();
    cY = tracks[i]->helix().center().y();
    if(fabs(r)*(1.-m_param.MERGE_RATIO) <= R && R <= fabs(r)*(1.+m_param.MERGE_RATIO)){
      if(fabs(cx-cX) <= fabs(r)*m_param.MERGE_RATIO && fabs(cy-cY) <= fabs(r)*m_param.MERGE_RATIO){
    double diff = fabs((fabs(r)-fabs(R))*(cx-cX)*(cy-cY));
    if(diff < bestDiff){
      bestDiff = diff;
      bestIndex = i;
    }
      }
    }
  }
  if(bestDiff == 1.0e20)return track;
  R  = tracks[bestIndex]->helix().radius();
  cX = tracks[bestIndex]->helix().center().x();
  cY = tracks[bestIndex]->helix().center().y();
  if(r*R >= 0.){
    if(fabs(track->helix().dz()-tracks[bestIndex]->helix().dz()) < m_param.MERGE_Z_DIFF){
      if(track->nLinks() > tracks[bestIndex]->nLinks()){
    m_tracks.remove(tracks[bestIndex]);
    return track;
      }else{
    m_tracks.remove(track);
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Merged...removing this track.(type1)" << std::endl;
#endif
    return NULL;
      }
    }
  }else{
    bool newTrack(false), oldTrack(false);
    unsigned newCounter(0), oldCounter(0);
    for(unsigned i=0, size=m_hitsOnInnerSuperLayer.length();
    i<size;++i){
      if(!oldTrack){
    if((R > 0. && 
        cX*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY)-
        cY*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX) > 0.) ||
       (R < 0. && 
        cX*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY)-
        cY*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX) < 0.)){
      double dist = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cX,
                 m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cY);
      if(dist < fabs(R)){
        if(fabs(fabs(R)-dist-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
          ++oldCounter;
          if(oldCounter >= 3)oldTrack = true;
        }
      }else{
        if(fabs(dist-fabs(R)-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
          ++oldCounter;
          if(oldCounter >= 3)oldTrack = true;
        }
      }
    }
      }
      if(!newTrack){
    if((r > 0. && 
        cx*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy)-
        cy*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx) > 0.) ||
       (r < 0. && 
        cx*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy)-
        cy*(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx) < 0.)){
      double dist = distance(m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().x()-cx,
                 m_hitsOnInnerSuperLayer[i]->wire()->xyPosition().y()-cy);
      if(dist < fabs(r)){
        if(fabs(fabs(r)-dist-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
          ++newCounter;
          if(newCounter >= 3)newTrack = true;
        }
      }else{
        if(fabs(dist-fabs(r)-m_hitsOnInnerSuperLayer[i]->drift()) < 0.5){
          ++newCounter;
          if(newCounter >= 3)newTrack = true;
        }
      }
    }
      }
      if(oldTrack && newTrack)break;
    }
    if(oldTrack && !newTrack){
      m_tracks.remove(track);
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)       3D:Merged...removing this track.(type2)" << std::endl;
#endif
      return NULL;
    }else if(!oldTrack && newTrack){
      m_tracks.remove(tracks[bestIndex]);
      return track;
    }else if(!oldTrack && !newTrack){
      m_tracks.remove(track);
      m_tracks.remove(tracks[bestIndex]);
#if DEBUG_CURL_DUMP
      std::cout << "(TCurlFinder)       3D:Merged...removing this track.(type3)" << std::endl;
#endif
      return NULL;
    }else if(oldTrack && newTrack){
      if(fabs(track->helix().dz()) > fabs(tracks[bestIndex]->helix().dz()) &&
     fabs(track->helix().dz()) > fabs(tracks[bestIndex]->helix().dz())+m_param.MERGE_Z_DIFF){
    m_tracks.remove(track);
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Merged...removing this track.(type4)" << std::endl;
#endif
    return NULL;
      }else if(fabs(tracks[bestIndex]->helix().dz()) > fabs(track->helix().dz()) &&
           fabs(tracks[bestIndex]->helix().dz()) > fabs(track->helix().dz())+m_param.MERGE_Z_DIFF){
    m_tracks.remove(tracks[bestIndex]);
    return track;
      }
    }
  }
  return track;
}
 
void
TCurlFinder::salvage3DTrack(TTrack *track, bool half) {
  if(track->nLinks() >= m_param.MIN_SALVAGE){
    AList<TMLink> list = m_unusedAxialHits;
    list.append(m_unusedStereoHits);
    list.remove(track->links());
 
    if(half){
      double q = track->charge();
      double x = track->helix().center().x();
      double y = track->helix().center().y();
      AList<TMLink> removeList;
      const double Bz = -1000*m_pmgnIMF->getReferField();
      for(unsigned i = 0, size = list.length(); i < size; ++i){
    if(q*Bz*(x*list[i]->wire()->xyPosition().y()-y*list[i]->wire()->xyPosition().x())<0.){
    //if(q*(x*list[i]->wire()->xyPosition().y()-y*list[i]->wire()->xyPosition().x())<0.)
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" removeList append  "<<i<<std::endl;
#endif
      removeList.append(list[i]);
    }
      }
      list.remove(removeList);
    }
 
    AList<TMLink> badCand, goodCand;
    double dist;
    for(unsigned j=0, nLinks=track->nLinks();j<nLinks;++j){
      dist = distance(*track, *(track->links()[j]));
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" "<<j<<" distance  "<<dist<<std::endl;
#endif
      if(dist > m_param.BAD_DISTANCE_FOR_SALVAGE)badCand.append(track->links()[j]);
    }
    for(unsigned j=0, nList=list.length();j<nList;++j){
      dist = distance(*track, *(list[j]));
#ifdef TRKRECO_DEBUG_DETAIL
      std::cout<<"   "<<__FILE__<<"   "<<__LINE__<<" "<<j<<" distance  "<<dist<<std::endl;
#endif
      if(dist < m_param.GOOD_DISTANCE_FOR_SALVAGE)goodCand.append(list[j]);
    }
    track->TTrackBase::remove(badCand);
    track->TTrackBase::append(goodCand);
    if(m_fitter.fit(*track) < 0){
      track->TTrackBase::remove(goodCand);
      track->TTrackBase::append(badCand);
      m_fitter.fit(*track);
    }
  }
  return;
}
 
double 
TCurlFinder::distance(const TTrack &track, const TMLink &link) const {
//  if(link.wire()->axial()){      Liuqg 060926
  if((link.wire()->superLayerId() > 1 && link.wire()->superLayerId() <5) || link.wire()->superLayerId() >8 ){
    //...axial
    double d = distance(track.helix().center().x()-link.xyPosition().x(),
            track.helix().center().y()-link.xyPosition().y());
    double diff = fabs(d - fabs(track.helix().radius()));
    return fabs(link.hit()->drift()-diff);
  }
  //...stereo
  HepPoint3D xc(track.helix().center());
  HepPoint3D xw(link.xyPosition());
  HepPoint3D xt(track.helix().x());
  HepVector3D v0(xt-xc), v1(xw-xc);
  double vCrs(v0.x() * v1.y() - v0.y() * v1.x());
  double vDot(v0.x() * v1.x() + v0.y() * v1.y());
  double dPhi = atan2(vCrs, vDot);
  Vector a(track.helix().a());
  double kappa = a[2];
  double phi0  = a[1];
 
  const double Bz = -1000*m_pmgnIMF->getReferField();
  //yzhang 2012-05-03 
  double rho = m_param.ALPHA_SAME_WITH_HELIX/kappa/Bz;
  double tanLambda = a[4];
  HepVector3D v = link.wire()->direction();
  Vector c(3);
  c = HepPoint3D(link.wire()->backwardPosition()-(v*link.wire()->backwardPosition())*v);
  
  HepDiagMatrix e(3,1);
  Matrix t(3, 3);
  t[0][0] = v.x() * v.x();
  t[0][1] = v.x() * v.y();
  t[0][2] = v.x() * v.z();
  t[1][0] = t[0][1];
  t[1][1] = v.y() * v.y();
  t[1][2] = v.y() * v.z();
  t[2][0] = t[0][2];
  t[2][1] = t[1][2];
  t[2][2] = v.z() * v.z();
  t -= e;
 
  double factor = 1.;
  unsigned nTrial = 0;
  
  //...Cal. closest point(Newton method)...
  Vector x(3);
  Vector dXdPhi(3);
  Vector d2Xd2Phi(3);
  double fOld = 0.; // The initialization is not needed.
  const double convergence = 1.0e-5;
  while(nTrial < 100){
    x = track.helix().x(dPhi);
    double cosPhi = cos(phi0+dPhi);
    double sinPhi = sin(phi0+dPhi);
    dXdPhi[0] =   rho*sinPhi;
    dXdPhi[1] = - rho*cosPhi;
    dXdPhi[2] = - rho*tanLambda;
    
    //...f = d(Distance) / d phi...
    double f = dot(c,dXdPhi)+dot(x,(t*dXdPhi));
    
    if(fabs(f) < convergence)break;
    if(nTrial > 0){
      double eval = (1.-0.25*factor)*fabs(fOld)-fabs(f);
      if(eval <= 0.)factor *= 0.5;
    }
    //...Cal. next phi...
    d2Xd2Phi[0] = rho*cosPhi;
    d2Xd2Phi[1] = rho*sinPhi;
    d2Xd2Phi[2] = 0.;
    double df = dot(c, d2Xd2Phi)+ 
      dot(dXdPhi, (t * dXdPhi))+
      dot(x, (t * d2Xd2Phi));
    dPhi -= factor * f / df;
    
    fOld = f;
    ++nTrial;
  }
 
  double beta = v*(track.helix().x(dPhi)-link.wire()->backwardPosition());
  return fabs((link.wire()->backwardPosition()+beta*v-track.helix().x(dPhi)).mag()- 
          link.hit()->drift());
}
 
//............................................
//.................2D Parts...................
//............................................
 
TCircle *
TCurlFinder::make2DTrack(const AList<TMLink> &seed,
             const AList<TSegmentCurl> &segmentList,
             const unsigned ip) {
  //jialk origin is 3            
  if(seed.length() < m_param.minimum_seedLength)return NULL;
  TCircle *circle = new TCircle(seed);
  m_allCircles.append(circle);
  int errorFlag = circle->fitForCurl(ip);
  //jialk add a limitation of ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK
  //if(fabs(circle->radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK)return NULL;
  if( (fabs(circle->radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK) || (fabs(circle->radius()) < m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK) )return NULL;
  int searchDirection = 1;//[ 1 : inner ], [ -1 : outer ]
  int searchPath(searchDirection);
  bool searchZero(false);
  bool changeDirection(false);
  unsigned superLayerId = seed[0]->hit()->wire()->superLayerId();
 
  AList<TMLink> cand, tmpList;
  AList<TMLink> preAxialCand, preStereoCand;
  makeList(tmpList, segmentList, seed);
  for(unsigned i=0, size=tmpList.length();i<size;++i){
    if(tmpList[i]->wire()){
      if(tmpList[i]->wire()->axial()) preAxialCand.append(tmpList[i]);
      else preStereoCand.append(tmpList[i]);
    }
  }
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)  2D: Superlayer of seed = " << superLayerId << std::endl;
#endif
 
  bool appendFlag = false;
nextStep:
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)  2D: SearchPath = " << searchPath
       << " Search SelfSuperlayer = " << (int)(searchZero)
       << " Change Direction of Search = " << (int)(changeDirection) << std::endl;
#endif
  if(preAxialCand.length() == 0 && preStereoCand.length() == 0){
    if(circle->links().length() >= 3){
      if(m_unusedAxialHits.length() == 0)
    if(errorFlag == -1)return NULL;
    else return circle;
      else goto salvage;
    }else{
      if(m_unusedAxialHits.length() == 0)return NULL;
      else goto salvage;      
    }
  }
  searchAxialCand(cand, preAxialCand, circle, 
          searchPath, superLayerId, m_param.RANGE_FOR_AXIAL_SEARCH);
  if(cand.length() > 0){
    appendFlag = true;
    for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
    errorFlag = circle->fitForCurl(ip);
    preAxialCand.remove(circle->links());
    if(m_param.STEREO_2DFIND){
      searchStereoCand(cand, preStereoCand, circle, 
               searchPath, superLayerId, m_param.RANGE_FOR_STEREO_SEARCH);
      if(cand.length() > 0){
    appendFlag = true;
    for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
    errorFlag = circle->fitForCurl(ip);
    preStereoCand.remove(circle->links());
      }
    }
    if(searchDirection == 1)++searchPath;
    else --searchPath;
    goto nextStep;
  }else{
    if(m_param.STEREO_2DFIND){
      searchStereoCand(cand, preStereoCand, circle, 
               searchPath, superLayerId, m_param.RANGE_FOR_STEREO_SEARCH);
      if(cand.length() > 0){
    appendFlag = true;
    for(unsigned i = 0, size = cand.length(); i < size; ++i)circle->append(*cand[i]);
    errorFlag = circle->fitForCurl(ip);
    preStereoCand.remove(circle->links());
    if(searchDirection == 1)++searchPath;
    else --searchPath;
    goto nextStep;
      }else if((searchPath == 1 || searchPath == -1) && !searchZero){
    searchPath = 0;
    searchZero = true;
    goto nextStep;
      }else if((searchPath == 1 || searchPath == -1) && searchZero && !changeDirection){
    searchPath *= -1;
    searchDirection *= -1;
    changeDirection = true;
    goto nextStep;
      }else{
    if(circle->links().length() >= 3){
      if(m_unusedAxialHits.length() == 0)
        if(errorFlag == -1)return NULL;
        else return circle;
      else goto salvage;
    }else{
      if(m_unusedAxialHits.length() == 0)return NULL;
      else goto salvage;
    }
      }
    }else{
      if((searchPath == 1 || searchPath == -1) && !searchZero){
    searchPath = 0;
    searchZero = true;
    goto nextStep;
      }else if((searchPath == 1 || searchPath == -1) && searchZero && !changeDirection){
    searchPath *= -1;
    searchDirection *= -1;
    changeDirection = true;
    goto nextStep;
      }else{
    if(circle->links().length() >= 3){
      if(m_unusedAxialHits.length() == 0)
        if(errorFlag == -1)return NULL;
        else return circle;
      else goto salvage;
    }else{
      if(m_unusedAxialHits.length() == 0)return NULL;
      else goto salvage;
    }
      }
    }
  }
 
salvage:
  cand.removeAll();
  searchHits(cand, m_unusedAxialHits,  circle, m_param.RANGE_FOR_AXIAL_LAST2D_SEARCH);
  if(m_param.STEREO_2DFIND){
    searchHits(cand, m_unusedStereoHits, circle, m_param.RANGE_FOR_STEREO_LAST2D_SEARCH);
  }
  if(checkAppendHits(circle->links(), cand)){
    circle->append(cand);
    if(circle->nLinks() >= 3)
      if(circle->fitForCurl(ip) == -1)return NULL;
      else return circle;
    else return NULL;
  }else if(circle->nLinks() >= 3){
    return circle;
  }else return NULL;
}
 
TMLink *
findIsolatedCloseHits(TMLink *link);
 
void
TCurlFinder::searchAxialCand(AList<TMLink> &cand,
                 const AList<TMLink> &preCand,
                 const TCircle *circle,
                 const int depth,
                 const unsigned superLayerID,
                 const double searchError) {
  cand.removeAll();
  int innerSuperLayerId = nextSuperAxialLayerId(superLayerID, depth);
  if(innerSuperLayerId < 0)return;
  for(unsigned i = 0, size = preCand.length(); i < size; ++i){
    if(preCand[i]->hit()->wire()->superLayerId() == 
       (static_cast<unsigned>(innerSuperLayerId))){
#if 0
      if(searchHits(preCand[i], circle, searchError))cand.append(preCand[i]);
#else
      if(searchHits(preCand[i], circle, searchError)){
    cand.remove(preCand[i]);
    cand.append(preCand[i]);
    TMLink * cand2 = findIsolatedCloseHits(preCand[i]);
    if(cand2){
      for(unsigned j = 0; j < size; ++j){
        if(preCand[j]->wire()->id() == cand2->wire()->id()){
          cand.remove(cand2);
          cand.append(cand2);
#if 0
          std::cout << "Axial Appending....";
          std::cout << " layerID = " << cand2->wire()->layerId();
          std::cout << " localID = " << cand2->wire()->localId() << std::endl;
          if(searchHits(cand2, circle, searchError)){
        std::cout << " But this can be added by default!" << std::endl;
          }else{
        std::cout << " Good!! this cannot be added by default!" << std::endl;
          }
#endif      
          break;
        }
      }
    }
      }
#endif
    }
  }
}
 
void
TCurlFinder::searchStereoCand(AList<TMLink> &cand,
                  const AList<TMLink> &preCand,
                  const TCircle *circle,
                  const int depth,
                  const unsigned superLayerID,
                  const double searchError) {
  cand.removeAll();
  int innerSuperLayerId = nextSuperStereoLayerId(superLayerID, depth);
  if(innerSuperLayerId < 0 || innerSuperLayerId > m_param.SUPERLAYER_FOR_STEREO_SEARCH)return;
  for(unsigned i = 0, size = preCand.length(); i < size; ++i){
    if(preCand[i]->hit()->wire()->superLayerId() == 
       (static_cast<unsigned>(innerSuperLayerId))){
      if(searchHits(preCand[i], circle, searchError))cand.append(preCand[i]);
    }
  }
}
 
unsigned
TCurlFinder::searchHits(const TMLink *link, const TCircle *circle, 
            const double searchError) const {
  // ...checks whether "link" can be added to circle.
  // ..."searchError" is length for checking.
  // ...returns 0 = error
  //            1 = no error
  double dist = distance(link->hit()->wire()->xyPosition().x() - circle->center().x(),
             link->hit()->wire()->xyPosition().y() - circle->center().y());
  double radius = fabs(circle->radius());
  //std::cout << link->wire()->localId() << " " << link->wire()->layerId() << " r=" 
  //   << radius << " e=" << searchError << " d=" << dist << std::endl;
  if(radius - searchError < dist &&
     radius + searchError > dist)return 1;
  return 0;
}
 
unsigned
TCurlFinder::searchHits(AList<TMLink> &cand,
            const AList<TMLink> &preCand,
            const TCircle *circle,
            const double searchError) const {
  unsigned numBefore = cand.length();
  for(unsigned i = 0, size = preCand.length(); i < size; ++i){
    if(searchHits(preCand[i], circle, searchError)){
#if 0
      cand.append(preCand[i]);
#else
      cand.remove(preCand[i]);
      cand.append(preCand[i]);
      TMLink * cand2 = findIsolatedCloseHits(preCand[i]);
      if(cand2){
    for(unsigned j = 0; j < size; ++j){
      if(preCand[j]->wire()->id() == cand2->wire()->id()){
        cand.remove(cand2);
        cand.append(cand2);
        break;
      }
    }
      }
#endif
    }
  }
  if(numBefore == cand.length())return 0;
  return 1;
}
 
unsigned
TCurlFinder::checkAppendHits(const AList<TMLink> &link, 
                 AList<TMLink> &cand) const {
  if(cand.length() == 0)return 0;
  AList<TMLink> tmp;
  for(unsigned i = 0, size1 = cand.length(), size2 = link.length(); i < size1; ++i){
    for(unsigned j = 0; j < size2; ++j){
      if((cand[i])->wire()->id() == (link[j])->wire()->id()){
    tmp.append(cand[i]);
    break;
      }
    }
  }
  cand.remove(tmp);
  if(cand.length() > 0)return 1;
  return 0;
}
 
void
TCurlFinder::removeStereo(TCircle &c) const
{
  AList<TMLink> stereoList;
  for(int i=0;i<c.links().length();++i){
    if(c.links()[i]->wire()->stereo()){
      stereoList.append(c.links()[i]);
    }
  }
  if(stereoList.length()>0)c.remove(stereoList);
}
 
bool
TCurlFinder::fitWDD(TCircle &c,
            double &chi2,
            int &ndf) const
{
  if(c.links().length() <= 3)return false;
  Lpav circle;
  // MDC
  for(int i=0;i<c.links().length();++i){
    circle.add_point((c.links()[i])->wire()->xyPosition().x(),     
                     (c.links()[i])->wire()->xyPosition().y(),1.0);
  }
  circle.add_point(0.,0.,1.0); // IP Constraint
  if (circle.fit() < 0.0 || circle.kappa() == 0.0) return false;
  double xc = circle.center()[0];
  double yc = circle.center()[1];
  double r  = circle.radius();
  const int maxIte = 2;
  for(int ite=0;ite<maxIte;++ite){
    Lpav circle2;
    circle2.clear();
    // MDC
    for(int i=0;i<c.links().length();++i){
      if(!((c.links()[i])->hit()->state() & WireHitFittingValid))continue;
      double R = sqrt(((c.links()[i])->wire()->xyPosition().x()-xc)*((c.links()[i])->wire()->xyPosition().x()-xc)+
                      ((c.links()[i])->wire()->xyPosition().y()-yc)*((c.links()[i])->wire()->xyPosition().y()-yc));
      if(R == 0.)continue;
      double U = 1./R;
      double dir = R > r ? -1. : 1.;
      double X = xc+((c.links()[i])->wire()->xyPosition().x()-xc)*U*(R+dir*(c.links()[i])->hit()->drift());
      double Y = yc+((c.links()[i])->wire()->xyPosition().y()-yc)*U*(R+dir*(c.links()[i])->hit()->drift());
      circle2.add_point(X,Y,1.0);
    }
    circle2.add_point(0.,0.,1.0); // IP Constraint
    if (circle2.fit() < 0.0 || circle2.kappa() == 0.0) return false;
    xc = circle2.center()[0];
    yc = circle2.center()[1];
    r  = circle2.radius();    
    //cout << xc << ", " << yc << " : " << r << endl;
  }
 
  // update of point information
  double totalChi2 = 0.;
  int totalNHit = 0;
  for(int i=0;i<c.links().length();++i){
    if(!((c.links()[i])->hit()->state() & WireHitFittingValid))continue;
    double xw = (c.links()[i])->wire()->xyPosition().x();
    double yw = (c.links()[i])->wire()->xyPosition().y();
    double R = sqrt((xw-xc)*(xw-xc)+(yw-yc)*(yw-yc));
    if(R == 0.)continue;
    double U = 1./R;
    double X = xc+(xw-xc)*U*r;
    double Y = yc+(yw-yc)*U*r;
    double zlr = xw*Y-yw*X;
    unsigned leftRight = zlr > 0. ? WireHitRight : WireHitLeft;
    double pChi2 = sqrt((X-xw)*(X-xw)+(Y-yw)*(Y-yw))-(c.links()[i])->hit()->drift();
    //cout << sqrt((X-xw)*(X-xw)+(Y-yw)*(Y-yw)) << " - " << (c.links()[i])->hit()->drift() << endl;
    //cout << i << ": " << pChi2 << endl;
    if((c.links()[i])->hit()->dDrift() != 0.){
      pChi2 *= pChi2/((c.links()[i])->hit()->dDrift()*(c.links()[i])->hit()->dDrift());
      totalChi2 += pChi2;
      //cout << pChi2 << ", " << c.links()[i]->hit()->dDrift() << endl;
      ++totalNHit;
    }else pChi2 = 1.0e+10;
    (c.links()[i])->update(HepPoint3D(X,Y,0.),HepPoint3D(xw,yw,0.),leftRight,pChi2);
    //cout << i << ": trk(" << X << "," << Y << "), wir(" << xw << "," << yw << ")" << endl;
  }
  chi2 = totalChi2;
  if(totalNHit <= 3)return false;
  ndf = totalNHit-3;
 
  HepPoint3D center(xc,yc,0.);
  double charge = 0.;
  //...Determine charge...Better way???
  int qSum = 0;
  for(int i=0;i<c.links().length();++i){
    TMLink * l = c.links()[i];
    if(l == 0)continue;   
    const TMDCWireHit * h = l->hit();
    if(h == 0)continue;    
    double q = (center.cross(h->xyPosition())).z();
    if(q > 0.)qSum += 1;
    else      qSum -= 1;
  }
  if(qSum >= 0)charge = +1.;
  else         charge = -1.;
  r *= charge;
  //cout << "B q = " << c.charge() << ", r = " << c.radius() << ", center = " << c.center() << endl;
  c.property(charge,r,center);
  //cout << "A q = " << c.charge() << ", r = " << c.radius() << ", center = " << c.center() << endl;
  return true;
}
 
//............................................
//.................3D Parts...................
//............................................
 
TTrack *
TCurlFinder::make3DTrack(const TCircle *circle, AList<TSegmentCurl> &segmentList) {
  unsigned size = segmentList.length();
  if(TTrack *track = make3DTrack(circle)){
    m_tracks.append(track);
#if 0
    std::cout << "MDC Helix+Pt: " << track->helix().dr() << ", " 
     << track->helix().phi0() << ", " 
     << track->helix().kappa() << ", " 
     << track->helix().dz() << ", " 
     << track->helix().tanl() 
     << ": " << 10000./2.9979258/15./track->helix().kappa() << std::endl;
#endif
    return track;
  }
  return NULL;
}
 
TTrack *
TCurlFinder::make3DTrack(const TCircle *circle) {
  TTrack *track = new TTrack(*circle);
  m_allTracks.append(track);
  //cout<<"2D: nAHits: "<<track->links().length()<<endl;
  //jialk caution origin is 3
  if(track->links().length() < m_param.minimum_2DTrackLength){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Fail...inital hit wire # < 3." << std::endl;
#endif
    return NULL;
  }
  AList<TMLink> allStereoHits(m_unusedStereoHits);
  allStereoHits.remove(track->links());
  AList<TMLink> closeHits;
  findCloseHits(allStereoHits, *track, closeHits);
  //if(!m_builder.buildStereo(*track, closeHits)){
  //jialk inorder to improve track robustness
    if ( closeHits.length() < m_param.minimum_closeHitsLength ){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Fail...stereohit wire # < 5." << std::endl;
#endif
    return NULL;
    }
  //if(!m_builder.buildStereo(*track, closeHits)){                    
  if(!m_builder.buildStereo(*track, closeHits, m_allStereoHitsOriginal)){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Fail...can not build stereo." << std::endl;
#endif
    return NULL;
  }
  //jialk add a limitation ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK 
  //if(fabs(track->helix().radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK){
  if( (fabs(circle->radius()) > m_param.MIN_RADIUS_OF_STRANGE_TRACK) || (fabs(circle->radius()) < m_param.ULTIMATE_MIN_RADIUS_OF_STRANGE_TRACK) ){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Fail...success 3D fit, but large radius > " 
     << m_param.MIN_RADIUS_OF_STRANGE_TRACK << "." << std::endl;
#endif
    return NULL;
  }
  //jialk caution origin is 5
  if(track->links().length() >= m_param.minimum_3DTrackLength){
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Success...can build stereo!!!" << std::endl;
#endif   
    return track;
  }else{
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)       3D:Fail...success 3D fit, but final hit wire # < 3." << std::endl;
#endif   
    return NULL;
  }
}
 
TMLink *
findIsolatedCloseHits(TMLink *link)
{
  int nNeighbor = 0;
  int nIsolated = 0;
  TMLink * isolatedLink[2] = {NULL,NULL};
  unsigned layerID = link->wire()->layerId();
  unsigned localID = link->wire()->localId();
  for(int i=0;i<6;++i){
    if(link->neighbor(i)){
      if(link->neighbor(i)->wire()->layerId() == layerID){
    ++nNeighbor;
    int isolated = 1;
    int testEach = 0; // for test
    for(int j=0;j<6;++j){
      if(link->neighbor(i)->neighbor(j)){
        if(link->neighbor(i)->neighbor(j)->wire()->layerId() == layerID &&
           link->neighbor(i)->neighbor(j)->wire()->localId() != localID){
          isolated = 0;
        }
#if 1
        else if(link->neighbor(i)->neighbor(j)->wire()->layerId() == layerID &&
            link->neighbor(i)->neighbor(j)->wire()->localId() == localID){
          testEach = 1;
        }
#endif
      }else break;
    }
    if(isolated == 1){
      if(nIsolated < 2)
        isolatedLink[nIsolated] = link->neighbor(i);
      ++nIsolated;    
    }
#if 1
    if(testEach == 0){
      std::cout << "Why?? Neighborhood info. dose not exist!!" << std::endl;
    }
#endif
      }
    }else break;
  }
#if 0
  std::cout << "isolated/neighbor # = " << nIsolated << "/" << nNeighbor << std::endl;
  std::cout << "layer ID = " << layerID << " ";
  std::cout << "local ID = " << localID << " --> ";
  if(isolatedLink[0])std::cout << isolatedLink[0]->wire()->localId() << " ";
  if(isolatedLink[1])std::cout << isolatedLink[1]->wire()->localId() << " ";
  std::cout << std::endl;
#endif
  if(nIsolated == 1 &&
     nNeighbor == 1 && isolatedLink[0])return isolatedLink[0];
  else return NULL;
}
 
void
TCurlFinder::findCloseHits(AList<TMLink> &links,
               TTrack &track, AList<TMLink> &list) {
  // ...finds candidates in the "links".
  // ...Candidates mean |"track" - wire(from "links" elements)| < dRcut 
  // ...returns these candidates("list").
  double dRcut[11] = {m_param.RANGE_FOR_STEREO_FIRST, m_param.RANGE_FOR_STEREO_SECOND, 
              0., 0., 
                      0., m_param.RANGE_FOR_STEREO_THIRD,
              m_param.RANGE_FOR_STEREO_FORTH, m_param.RANGE_FOR_STEREO_FIFTH,
              m_param.RANGE_FOR_STEREO_SIXTH, 0.,
                      0.};
  double r = fabs(track.helix().curv());
  double q = track.charge();
  double x = track.helix().center().x();
  double y = track.helix().center().y();
  for(unsigned i = 0, size = links.length(); i < size; ++i){
    if(fabs((links[i]->wire()->xyPosition() - track.helix().center()).mag() - r) < 
       dRcut[links[i]->wire()->superLayerId()]){
      if(q*(x*links[i]->wire()->xyPosition().y()-y*links[i]->wire()->xyPosition().x()) > 0.){
    list.remove(links[i]);
    list.append(links[i]);
    TMLink *cand = findIsolatedCloseHits(links[i]);
    if(cand){
      list.remove(cand);
      list.append(cand);
    }
      }
    }
  }
  return;
}
 
//..................................................
//..................................................
//..................................................
//..................................................
//..................................................
 
int 
TCurlFinder::makeWithMC(const AList<TMDCWireHit> & axialHits,
            const AList<TMDCWireHit> & stereoHits,
            AList<TTrack> & tracks) {
#if DEBUG_CURL_MC
#define MAX_INDEX_MAKEMC 100
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)Now making tracks using MC info..." << std::endl;
#endif
  int index[MAX_INDEX_MAKEMC];
  for(unsigned i = 0; i < MAX_INDEX_MAKEMC; ++i)index[i] = 9999;
 
  int counter(0);
  bool first(true);
 
  for(unsigned i = 0, size = axialHits.length(); i < size; ++i){
    if(axialHits[i]->mc() &&
       axialHits[i]->mc()->hep() &&
       !(axialHits[i]->state() & WireHitUsed)){
      int flag(1);
      for(unsigned j = 0; j < MAX_INDEX_MAKEMC; ++j){
    if(index[j] != 9999 && index[j] == axialHits[i]->mc()->hep()->id()){
      flag = 0;
      break;
    }
      }
      if(flag){
    index[counter] = axialHits[i]->mc()->hep()->id();
    ++counter;
      }
    }
  }
#if DEBUG_CURL_DUMP
  std::cout << "(TCurlFinder)Found " << counter 
       << " tracks with MC information." << std::endl;
#endif
  for(unsigned j = 0; j < counter; ++j){
    AList<TMLink> axialList;
    AList<TMLink> stereoList;
    int axialCounter(0);
    int stereoCounter(0);
    //...axial
    for(unsigned i = 0, size = axialHits.length(); i < size; ++i){
      if(index[j] == axialHits[i]->mc()->hep()->id() &&
     !(axialHits[i]->state() & WireHitUsed)){
    axialList.append(new TMLink(0, axialHits[i]));
    ++axialCounter;
      }
    }
    if(axialCounter < 3){
      HepAListDeleteAll(axialList);
      continue;
    }
    //...stereo
    for(unsigned i = 0, size = stereoHits.length(); i < size; ++i){
      if(index[j] == stereoHits[i]->mc()->hep()->id() &&
     !(stereoHits[i]->state() & WireHitUsed)){
    stereoList.append(new TMLink(0, stereoHits[i]));
    ++stereoCounter;
      }
    }
    if(stereoCounter < 2){
      HepAListDeleteAll(axialList);
      HepAListDeleteAll(stereoList);
      continue;
    }
#if DEBUG_CURL_DUMP
    std::cout << "(TCurlFinder)#" << j << " : Use " 
     << axialCounter << " axial hit wires and "
     << stereoCounter << " stereo hit wires" << std::endl;
    std::cout << "(TCurlFinder)Particle Type(LUND) = " 
     << axialList[0]->hit()->mc()->hep()->pType() << std::endl;
#endif
 
    m_unusedAxialHitsOriginal.append(axialList);
    m_unusedAxialHitsOriginal.append(stereoList);
    TCircle *circle = new TCircle(axialList);
    m_allCircles.append(circle);
    circle->fitForCurl();
    double charge = 1.;
    if(axialList[0]->hit()->mc()->hep()->pType() < 0)charge = -1.;
    if(fabs(axialList[0]->hit()->mc()->hep()->pType()) == 11 ||
       fabs(axialList[0]->hit()->mc()->hep()->pType()) == 13 ||
       fabs(axialList[0]->hit()->mc()->hep()->pType()) == 15)charge *= -1.;
    circle->property(charge, charge*fabs(circle->radius()), circle->center());
 
    AList<TMLink> removeList;
    double x = circle->center().x();
    double y = circle->center().y();
    //...axial
    for(unsigned i = 0, size = axialList.length();
    i < size; ++i){
      if(charge*(x*axialList[i]->xyPosition().y()-
         y*axialList[i]->xyPosition().x())< 0.){
    removeList.append(axialList[i]);
 
      }
    }
    circle->remove(removeList);
    if(circle->nLinks() < 3)continue;
    //...refits
    circle->fitForCurl(1);
    x = circle->center().x();
    y = circle->center().y();
    removeList.removeAll();
    ///...stereo
    for(unsigned i = 0, size = stereoList.length();
    i < size; ++i){
      if(charge*(x*stereoList[i]->xyPosition().y()-
         y*stereoList[i]->xyPosition().x())< 0.){
    removeList.append(stereoList[i]);
      }
    }
    stereoList.remove(removeList);
    if(stereoList.length() < 2)continue;
 
    TTrack *track = new TTrack(*circle);
    m_allTracks.append(track);
    if(m_builder.buildStereoMC(*track, stereoList)){
      if(track->links().length() >= 5){
    track->assign(WireHitCurlFinder, TrackCurlFinder | TrackValid | Track3D);
    tracks.append(track);
    m_allTracks.remove(track);
      }else{
    std::cout << "Can not reconstruct with MC information!" << std::endl;
      }
    }else{
      std::cout << "Can not reconstruct with MC information!" << std::endl;
    }
  }
#endif
  return 0;
}
 
void
TCurlFinder::makeCdcFrame(void) {
#if DEBUG_CURL_GNUPLOT+DEBUG_CURL_SEGMENT
  //#if 1
  double X = 0.;
  double Y = 0.;
  double R[12] = {8.3,  16.9, 21.7, 31.3, 36.1, 44.1,
          50.5, 58.5, 64.9, 72.9, 79.3, 87.4};
  double step = 300.;
  double dStep = 2.*M_PI/step;
  FILE *data;
  std::string nameHead = "tmp.cdc_";
  for(int j=0;j<12;++j){
    std::string nameFile = nameHead+"0"+itostring(j);
    if(j>=10)nameFile = nameHead+itostring(j);
    if((data = fopen(nameFile,"w")) != NULL){
      for(int i=0;i<step;++i){
    double x = X + R[j] * cos(dStep*static_cast<double>(i));
    double y = Y + R[j] * sin(dStep*static_cast<double>(i));
    fprintf(data,"%lf, %lf\n",x,y);
      }
      fclose(data);
    }
  }
 
  if((data = fopen("tmp_wires.dat","w")) != NULL){
    AList<TMLink> list = m_unusedAxialHitsOriginal;
    list.append(m_unusedStereoHitsOriginal);
    for(int i=0;i<list.length();i++){
      double x = list[i]->hit()->wire()->xyPosition().x();
      double y = list[i]->hit()->wire()->xyPosition().y();
      fprintf(data,"%lf, %lf\n",x,y);
    }
    fclose(data);
  }
#endif
  return;
}
 
void 
TCurlFinder::plotSegment(const AList<TMLink>& list, const int flag) {
#if DEBUG_CURL_GNUPLOT
  if(!m_debugCdcFrame){
    makeCdcFrame();
    m_debugCdcFrame = true;
  }
  double gmaxX = 90. ,gminX = -90.;
  double gmaxY = 90. ,gminY = -90.;
  FILE *gnuplot, *data;
  if((data = fopen("tmp.dat","w")) != NULL){
    if(flag)std::cout << "Wire ID = ";
    for(int i=0;i<list.length();i++){
      double x = list[i]->hit()->wire()->xyPosition().x();
      double y = list[i]->hit()->wire()->xyPosition().y();
      fprintf(data,"%lf, %lf\n",x,y);
      if(flag)std::cout << list[i]->hit()->wire()->id() << ", ";
    }
    if(flag)std::cout << std::endl;
    fclose(data);
  }
  if((gnuplot = popen("gnuplot","w")) != NULL){
    fprintf(gnuplot,"set nokey \n");
    fprintf(gnuplot,"set size 0.721,1.0 \n");
    fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
    fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
    std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat\"";
    std::string nameHead = ",\"tmp.cdc_";
    for(int j=0;j<12;++j){
      std::string nameFile = nameHead+"0"+itostring(j)+"\"w l 0";
      if(j>=10)nameFile = nameHead+itostring(j)+"\"w l 0";
      longName += nameFile;
    }
    longName += " \n";
    fprintf(gnuplot,longName);
    fflush(gnuplot);
    char tmp[8];
    gets(tmp);
    pclose(gnuplot);
  }
#endif
  return;
}
 
void 
TCurlFinder::plotCircle(const TCircle& circle, const int flag) {
#if DEBUG_CURL_GNUPLOT
  //#if 1
  if(!m_debugCdcFrame){
    makeCdcFrame();
    m_debugCdcFrame = true;
  }
  double gmaxX = 90. ,gminX = -90.;
  double gmaxY = 90. ,gminY = -90.;
  FILE *gnuplot, *data;
  if((data = fopen("tmp.dat1","w")) != NULL){
    if(flag)std::cout << "Axial  Wire ID ==> " << std::endl;
    for(int i=0;i<circle.nLinks();++i){
      if(circle.links()[i]->hit()->wire()->axial()){
    double x = circle.links()[i]->hit()->wire()->xyPosition().x();
    double y = circle.links()[i]->hit()->wire()->xyPosition().y();
    fprintf(data,"%lf, %lf\n",x,y);
    if(flag){
      /*if(debugMcFlag){
        std::cout << " A:" << circle.links()[i]->hit()->wire()->id() << ", ";
        std::cout << ", HepTrackID = " << circle.links()[i]->hit()->mc()->hep()->id();
        std::cout << ", HepLundID = "  << circle.links()[i]->hit()->mc()->hep()->pType() << std::endl;
        }else std::cout << " A:" << circle.links()[i]->hit()->wire()->id() << std::endl;*/
    }
      }
    }
    if(flag)std::cout << std::endl;
    fclose(data);
  }
  if((data = fopen("tmp.dat2","w")) != NULL){
    if(flag)std::cout << "Stereo Wire ID ==> " << std::endl;
    for(int i=0;i<circle.nLinks();++i){
      if(circle.links()[i]->hit()->wire()->stereo()){
    double x = circle.links()[i]->hit()->wire()->xyPosition().x();
    double y = circle.links()[i]->hit()->wire()->xyPosition().y();
    fprintf(data,"%lf, %lf\n",x,y);
    if(flag){
      /*if(debugMcFlag){
        std::cout << " S:" << circle.links()[i]->hit()->wire()->id() << ", ";
        std::cout << ", HepTrackID = " << circle.links()[i]->hit()->mc()->hep()->id();
        std::cout << ", HepLundID = "  << circle.links()[i]->hit()->mc()->hep()->pType() << std::endl;
        }else std::cout << " S:" << circle.links()[i]->hit()->wire()->id() << std::endl;*/
    }
      }
    }
    if(flag)std::cout << std::endl;
    fclose(data);
  }
  double X = circle.center().x();
  double Y = circle.center().y();
  double R = fabs(circle.radius());
  double step = 300.;
  double dStep = 2.*M_PI/step;
  if((data = fopen("tmp.dat3","w")) != NULL){
    for(int i=0;i<step;++i){
      double x = X + R * cos(dStep*static_cast<double>(i));
      double y = Y + R * sin(dStep*static_cast<double>(i));
      fprintf(data,"%lf, %lf\n",x,y);
    }
    fclose(data);
  }
  if((gnuplot = popen("gnuplot","w")) != NULL){
    fprintf(gnuplot,"set nokey \n");
    fprintf(gnuplot,"set size 0.721,1.0 \n");
    fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
    fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
    std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat1\", \"tmp.dat3\" w l, \"tmp.dat2\"";
    std::string nameHead = ",\"tmp.cdc_";
    for(int j=0;j<12;++j){
      std::string nameFile = nameHead+"0"+std::string(j)+"\"w l 0";
      if(j>=10)nameFile = nameHead+std::string(j)+"\"w l 0";
      longName += nameFile;
    }
    longName += " \n";
    fprintf(gnuplot,longName);
    fflush(gnuplot);
    char tmp[8];
    gets(tmp);
    pclose(gnuplot);
  }
#endif
  return;
}
 
void 
TCurlFinder::plotTrack(const TTrack& track, const int flag) {
#if DEBUG_CURL_GNUPLOT
  if(!m_debugCdcFrame){
    makeCdcFrame();
    m_debugCdcFrame = true;
  }
  double gmaxX = 90. ,gminX = -90.;
  double gmaxY = 90. ,gminY = -90.;
  FILE *gnuplot, *data;
  if((data = fopen("tmp.dat1","w")) != NULL){
    if(flag)std::cout << "Axial  Wire ID ==> " << std::endl;
    for(int i=0;i<track.nLinks();++i){
      if(track.links()[i]->hit()->wire()->axial()){
    double x = track.links()[i]->hit()->wire()->xyPosition().x();
    double y = track.links()[i]->hit()->wire()->xyPosition().y();
    fprintf(data,"%lf, %lf\n",x,y);
    if(flag){
      if(debugMcFlag){
        std::cout << " A:" << track.links()[i]->hit()->wire()->id() << ", ";
        std::cout << ", HepTrackID = " << track.links()[i]->hit()->mc()->hep()->id();
        std::cout << ", HepLundID = "  << track.links()[i]->hit()->mc()->hep()->pType() << std::endl;
      }else std::cout << " A:" << track.links()[i]->hit()->wire()->id() << std::endl;
    }
      }
    }
    if(flag)std::cout << std::endl;
    fclose(data);
  }
  if((data = fopen("tmp.dat2","w")) != NULL){
    if(flag)std::cout << "Stereo Wire ID ==> " << std::endl;
    for(int i=0;i<track.nLinks();++i){
      if(track.links()[i]->hit()->wire()->stereo()){
    double x = track.links()[i]->hit()->wire()->xyPosition().x();
    double y = track.links()[i]->hit()->wire()->xyPosition().y();
    fprintf(data,"%lf, %lf\n",x,y);
    if(flag){
      if(debugMcFlag){
        std::cout << " S:" << track.links()[i]->hit()->wire()->id() << ", ";
        std::cout << ", HepTrackID = " << track.links()[i]->hit()->mc()->hep()->id();
        std::cout << ", HepLundID = "  << track.links()[i]->hit()->mc()->hep()->pType() << std::endl;
      }else std::cout << " S:" << track.links()[i]->hit()->wire()->id() << std::endl;
    }
      }
    }
    if(flag)std::cout << std::endl;
    fclose(data);
  }
  double X = track.helix().center().x();
  double Y = track.helix().center().y();
  double R = fabs(track.helix().radius());
  double step = 300.;
  double dStep = 2.*M_PI/step;
  if((data = fopen("tmp.dat3","w")) != NULL){
    for(int i=0;i<step;++i){
      double x = X + R * cos(dStep*static_cast<double>(i));
      double y = Y + R * sin(dStep*static_cast<double>(i));
      fprintf(data,"%lf, %lf\n",x,y);
    }
    fclose(data);
  }
  if((gnuplot = popen("gnuplot","w")) != NULL){
    fprintf(gnuplot,"set nokey \n");
    fprintf(gnuplot,"set size 0.721,1.0 \n");
    fprintf(gnuplot,"set xrange [%f:%f] \n",gminX,gmaxX);
    fprintf(gnuplot,"set yrange [%f:%f] \n",gminY,gmaxY);
    std::string longName = "plot \"tmp_wires.dat\", \"tmp.dat1\", \"tmp.dat3\" w l, \"tmp.dat2\"";
    std::string nameHead = ",\"tmp.cdc_";
    for(int j=0;j<12;++j){
      std::string nameFile = nameHead+"0"+itostring(j)+"\"w l 0";
      if(j>=10)nameFile = nameHead+itostring(j)+"\"w l 0";
      longName += nameFile;
    }
    longName += " \n";
    fprintf(gnuplot,longName);
    fflush(gnuplot);
    char tmp[8];
    gets(tmp);
    pclose(gnuplot);
  }
#endif
  return;
}
 
void 
TCurlFinder::writeSegment(const AList<TMLink>& list, const int type) {
#if DEBUG_CURL_SEGMENT
  if(!m_debugCdcFrame){
    makeCdcFrame();
    m_debugCdcFrame = true;
  }
  double gmaxX = 90. ,gminX = -90.;
  double gmaxY = 90. ,gminY = -90.;
 
  FILE *data;
  std::string nameHead = "tmp.segment_";
  std::string nameFile = nameHead+itostring(type)+"_"+itostring(m_debugFileNumber);
  ++m_debugFileNumber;
  if((data = fopen(nameFile,"w")) != NULL){
    for(int i=0;i<list.length();i++){
      double x = list[i]->hit()->wire()->xyPosition().x();
      double y = list[i]->hit()->wire()->xyPosition().y();
      fprintf(data,"%lf, %lf\n",x,y);
    }
    fclose(data);
  }
#endif
  return;
}
 
void
TCurlFinder::dumpType1(TTrack *track) {
#if DEBUG_CURL_DUMP
  for(int j=0;j<track->nLinks();++j){
    std::cout << "Used Wire Info...";
    if(track->links()[j]->hit()->wire()->axial()){
      std::cout << "A:" << track->links()[j]->hit()->wire()->id() << ", ";
    }else{
      std::cout << "S:" << track->links()[j]->hit()->wire()->id() << ", ";
    }
    if(debugMcFlag){
      std::cout << ", HepTrackID = " << track->links()[j]->hit()->mc()->hep()->id();
      std::cout << ", HepLundID = "  << track->links()[j]->hit()->mc()->hep()->pType();
    }
    double dist = distance(*track, *(track->links()[j]));
    if(dist > 2.)std::cout << ": Large Distance( >2cm ) = " << dist;
    std::cout << std::endl;
  }
  AList<TMLink> list=m_unusedAxialHits;
  list.append(m_unusedStereoHits);
  for(unsigned j=0, nList=list.length();j<nList;++j){
    double dist = distance(*track, *(list[j]));
    std::cout << "Close Wire Info in ALL( <0.5cm )...";
    if(dist < 0.5){
      if(list[j]->hit()->wire()->axial())
    std::cout << "CA:" << list[j]->hit()->wire()->id() << ", ";
      else
    std::cout << "CS:" << list[j]->hit()->wire()->id() << ", ";
      if(debugMcFlag){
    std::cout << ", HepTrackID = " << list[j]->hit()->mc()->hep()->id();
    std::cout << ", HepLundID = "  << list[j]->hit()->mc()->hep()->pType();
      }
      std::cout << ", Distance = " << dist << std::endl;
    }
  }
#endif
  return;
}
 
void
TCurlFinder::dumpType2(TTrack *track) {
#if DEBUG_CURL_DUMP
  unsigned size = track->nLinks();
  if(size == 0)return;
 
  set< int, less<int> > uniqueHepID;
  vector<int> hepID;
  vector<double> ratio;
  for(int i=0;i<size;++i){
    uniqueHepID.insert(track->links()[i]->hit()->mc()->hep()->id());
    hepID.push_back(track->links()[i]->hit()->mc()->hep()->id());
    //std::cout << i << " : " << track->links()[i]->hit()->mc()->hep()->id() << std::endl;
  }
 
  set< int, less<int> >::iterator u = uniqueHepID.begin();
  vector<int>::size_type sizeInt;
  for(unsigned i=0;i<uniqueHepID.size();++i){
    sizeInt = 0;
    count(hepID.begin(), hepID.end(), *u, sizeInt);
    ratio.push_back((static_cast<double>(sizeInt)/static_cast<double>(size)));
    //std::cout << "HepID = " << *u << ", Ratio = " << ratio[i] << " = " << sizeInt << "/" << size << std::endl;
    ++u;
  }
 
  vector<double>::iterator m = max_element(ratio.begin(), ratio.end());
  int maxIndex = 0;
  ::distance(ratio.begin(), m, maxIndex);
  u = uniqueHepID.begin();
  advance(u,maxIndex);
  //std::cout << "MAX HepID = " << *u << ", Ratio = " << ratio[maxIndex] << std::endl;
  std::cout << "Ratio " << ratio[maxIndex] << std::endl;
  for(int i=0;i<size;++i){
    if(track->links()[i]->hit()->wire()->axial())std::cout << "A ";
    else std::cout << "S ";
 
    double dist = distance(*track, *(track->links()[i]));
    if(*u != track->links()[i]->hit()->mc()->hep()->id()){
      std::cout << "Bad " << dist << std::endl;
    }else{
      std::cout << "Good " << dist << std::endl;
    }
  }
#endif
  return;
}
 
//#if DEBUG_CURL_MC
#if 0
// --> TSegmentUtil.h
/*
sortBySvdRLA( const Datsvd_hit **a, const Datsvd_hit **b ) {
  if( (*a)->rla() < (*b)->rla() ){
    return -1;
  }else if( (*a)->rla() == (*b)->rla() ){
    return 0;
  }else{
    return 1;
  }
}
*/
 
Gen_hepevt *
cluster2hep(Recsvd_cluster *clus) {
  AList<Datsvd_hit> m_datsvd_hit;
  Datsvd_hit_Manager& svdHitMgr = Datsvd_hit_Manager::get_manager();
  for(Datsvd_hit_Manager::iterator
        it  = svdHitMgr.begin(),
        end = svdHitMgr.end();
      it != end; ++it){
    m_datsvd_hit.append(it);
  }
//  m_datsvd_hit.sort(sortBySvdRLA);
 
  unsigned size = m_datsvd_hit.length();
  int startPosition = -1;
  int direction = 1;
#if 0
  for(unsigned i=0;i<size;++i){
    //if(clus->width() == 1)
    std::cout << i << ": " << clus->hit().get_ID() << " <--> " << m_datsvd_hit[i]->get_ID()
     << ", RLA = " << m_datsvd_hit[i]->rla() << ", LSA = " << m_datsvd_hit[i]->amp()
     << ", Width = " << clus->width()
     << ", Cluster LSA = " << clus->lsa() << std::endl;
  }
#endif
  for(unsigned i=0;i<size;++i){
    if(m_datsvd_hit[i]->amp() == 0)std::cout << "DatSVD_Hit:amp == 0" << std::endl;
    //if(m_datsvd_hit[i]->amp() == 640)std::cout << "DatSVD_Hit:amp == 640" << std::endl;
    if(m_datsvd_hit[i]->rla() == 0)std::cout << "DatSVD_Hit:rla == 0" << std::endl;
    if(m_datsvd_hit[i]->rla() == 81920)std::cout << "DatSVD_Hit:rla == 81920" << std::endl;
    if(clus->hit().get_ID() == m_datsvd_hit[i]->get_ID()){
      startPosition = i;
      if(static_cast<double>(m_datsvd_hit[i]->amp()-1) > clus->lsa())direction = -1;
      break;
    }
  }
  if(startPosition == -1)return NULL;
  int width = clus->width();
#if 0
  std::cout << "Start # = " << startPosition 
       << ", Width = " << width 
       << ", Direction = " << direction <<std::endl;
#endif
 
  int* hepID = new int[width];
  set< int, less<int> > uniqueHepID;
 
  for(int i=startPosition;i<startPosition+width;++i)hepID[i-startPosition] = 0;
  //Datsvd_mchit_Manager& svdMcHitMgr = Datsvd_mchit_Manager::get_manager();
  Datsvd_mcdata_Manager& svdMcHitMgr = Datsvd_mcdata_Manager::get_manager();
  //std::cout << "SVD MC# = " << svdMcHitMgr.count() << std::endl;
  if(direction == 1){
    for(int i=startPosition;i<startPosition+width;++i){
      for(Datsvd_mcdata_Manager::iterator
        it  = svdMcHitMgr.begin(),
        end = svdMcHitMgr.end();
      it != end; ++it){
    if(it->Hep()){
      if(m_datsvd_hit[i]->rla() == it->rla() &&
         it->Hep().get_ID() != 0){
        hepID[i-startPosition] = it->Hep().get_ID();
        uniqueHepID.insert(it->Hep().get_ID());
        break;
      }
    }
      }
    }
  }else{
    int reverse = 0;
    for(int i=startPosition;i<startPosition+width;++i){
      ++reverse;
      //std::cout << startPosition+1-reverse << " --- " 
      //   << i << std::endl;
      for(Datsvd_mcdata_Manager::iterator
        it  = svdMcHitMgr.begin(),
        end = svdMcHitMgr.end();
      it != end; ++it){
    if(it->Hep()){
      if(m_datsvd_hit[startPosition+1-reverse]->rla() == it->rla() &&
         it->Hep().get_ID() != 0){
        hepID[i-startPosition] = it->Hep().get_ID();
        uniqueHepID.insert(it->Hep().get_ID());
        break;
      }
    }
      }
    }
  }
 
  unsigned num = uniqueHepID.size();
  int* counter = new int[num];
  set< int, less<int> >::iterator u = uniqueHepID.begin();
  for(int i=0;i<num;++i) counter[i] = 0;
 
  for(int i=0;i<num;++i){
    for(int j=0;j<width;++j){
      if(*u == hepID[j]){
    counter[i] += 1;
      }
    }
    ++u;
  }
 
  Gen_hepevt_Manager& genMgr = Gen_hepevt_Manager::get_manager();
#if 0
  u = uniqueHepID.begin();
  for(int i=0;i<num;++i){
    std::cout << i << ": TrackID = "<< *u - 1 
     << ", Count = " << counter[i]
     << ", LundID = " << genMgr[*u-1].idhep() << std::endl;
    ++u;
  }
#endif
 
  delete [] hepID;
  delete [] counter;
  if(num == 1){
    //std::cout << *(uniqueHepID.begin())-1 << std::endl;
    return &(genMgr[*(uniqueHepID.begin())-1]);
  }else if(num >= 2){
    //std::cout << "A svd cluster is not unique." << std::endl;
    return NULL;
  }else{
    return NULL;
  }
}
#endif
 
void
TCurlFinder::debugCheckSegments1(void) {
#if DEBUG_CURL_SEGMENT
  // Slow Checker(CPU time increases!!)
  // Neighboring wires should be included in the same segement.
  std::cout << "(TCurlFinder)checking consistency of segement..." << std::endl;
  unsigned nA = m_unusedAxialHitsOriginal.length();
  if(nA >= 2){
    for(unsigned i=0;i<nA-1;++i){
      int superLayerId = (int)(m_unusedAxialHitsOriginal[i]->wire()->superLayerId());
      int layerId  = (int)(m_unusedAxialHitsOriginal[i]->wire()->layerId());
      int localId  = (int)(m_unusedAxialHitsOriginal[i]->wire()->localId());
      int localIdP = (int)(m_unusedAxialHitsOriginal[i]->wire()->localIdForPlus());
      int localIdM = (int)(m_unusedAxialHitsOriginal[i]->wire()->localIdForMinus());
      for(unsigned j=i+1;j<nA;++j){
    int superLayerId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->superLayerId());
    int layerId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->layerId());
    int localId2 = (int)(m_unusedAxialHitsOriginal[j]->wire()->localId());
    if(superLayerId == superLayerId2){
      if(layerId2 == layerId){
        if(localIdP+1 == localId2 || localIdM-1 == localId2)
          debugCheckSegments(localId, layerId,
                 localId2,layerId2);
      }else if(layerId2 == layerId-1 || layerId2 == layerId+1){
        if(offset(layerId) == offset(layerId2)){
          std::cout << "(TCurlFinder: Waring) Offset is same at the same superlayer!!" << std::endl;
        }else if(offset(layerId) > offset(layerId2)){
          if(localId == localId2 || localIdP+1 == localId2)
        debugCheckSegments(localId, layerId,
                   localId2,layerId2);
        }else{
          if(localId == localId2 || localIdM-1 == localId2)
        debugCheckSegments(localId, layerId,
                   localId2,layerId2);
        }
      }
    }
      }
    }
  }
  unsigned nS = m_unusedStereoHitsOriginal.length();
  if(nS >= 2){
    for(unsigned i=0;i<nS-1;++i){
      int superLayerId = (int)(m_unusedStereoHitsOriginal[i]->wire()->superLayerId());
      int layerId  = (int)(m_unusedStereoHitsOriginal[i]->wire()->layerId());
      int localId  = (int)(m_unusedStereoHitsOriginal[i]->wire()->localId());
      int localIdP = (int)(m_unusedStereoHitsOriginal[i]->wire()->localIdForPlus());
      int localIdM = (int)(m_unusedStereoHitsOriginal[i]->wire()->localIdForMinus());
      for(unsigned j=i+1;j<nS;++j){
    int superLayerId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->superLayerId());
    int layerId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->layerId());
    int localId2 = (int)(m_unusedStereoHitsOriginal[j]->wire()->localId());
    if(superLayerId == superLayerId2){
      if(layerId2 == layerId){
        if(localIdP+1 == localId2 || localIdM-1 == localId2)
          debugCheckSegments(localId, layerId,
                 localId2,layerId2);
      }else if(layerId2 == layerId-1 || layerId2 == layerId+1){
        if(offset(layerId) == offset(layerId2)){
          std::cout << "(TCurlFinder: Waring) Offset is same at the same superlayer!!" << std::endl;
        }else if(offset(layerId) > offset(layerId2)){
          if(localId == localId2 || localIdP+1 == localId2)
        debugCheckSegments(localId, layerId,
                   localId2,layerId2);
        }else{
          if(localId == localId2 || localIdM-1 == localId2)
        debugCheckSegments(localId, layerId,
                   localId2,layerId2);
        }
      }
    }
      }
    }
  }
  std::cout << "(TCurlFinder)...done check of segement!" << std::endl;
  std::cout << "(TCurlFinder)...If no warning message exists, check of segement is complete!" << std::endl;
  std::cout << "(TCurlFinder)...Note: a segment size should be 1 or 2 to use this debugger." << std::endl;
#endif
  return;
}
 
void
TCurlFinder::debugCheckSegments(const double localId, const double layerId,
                const double localId2,const double layerId2) {
#if DEBUG_CURL_DUMP
  unsigned nSeg = m_segmentList.length();
  unsigned nFound = 0;
  for(unsigned i=0;i<nSeg;++i){
    unsigned nWire = m_segmentList[i]->list().length();
    unsigned mFound = 0;
    for(unsigned j=0;j<nWire;++j){
      if(((m_segmentList[i]->list())[j])->wire()->layerId() == layerId &&
     ((m_segmentList[i]->list())[j])->wire()->localId() == localId)++mFound;
      if(((m_segmentList[i]->list())[j])->wire()->layerId() == layerId2 &&
     ((m_segmentList[i]->list())[j])->wire()->localId() == localId2)++mFound;
    }
    if(mFound != 0 && mFound != 2){
      std::cout << "(TCurlFinder: Warning) Segment is inconsistency(0)!! mFound = " << mFound << std::endl;
    }
    if(mFound == 2)++nFound;
  }
  if(nFound != 1)
    std::cout << "(TCurlFinder: Warning) Segment is inconsistency(1)!! nFound = " << nFound << std::endl;
#endif
  return;
}
 
void
TCurlFinder::debugCheckSegments0(void) {
#if DEBUG_CURL_SEGMENT
  unsigned nSeg = m_segmentList.length();
  unsigned nWire = 0;
  for(unsigned i=0;i<nSeg;++i)nWire += m_segmentList[i]->list().length();
 
  unsigned nWireOriginal = m_unusedAxialHitsOriginal.length()+
                           m_unusedStereoHitsOriginal.length();
 
  std::cout << "(TCurlFinder: SelfChecker) Segment Parts" << std::endl;
  std::cout << "                           MIN_SEGMENT = " << m_param.MIN_SEGMENT << std::endl;
  std::cout << "                           Wire # of Orinal List = " << nWireOriginal << std::endl;
  std::cout << "                           Wire # of Segments    = " << nWire << std::endl;
  std::cout << "                           If MIN_SEGMENT <= 1, above numbers should be same." << std::endl;
  std::cout << "                           If MIN_SEGMENT >  1, former >= latter." << std::endl;
#endif
  return;
}
 
void
TCurlFinder::debugCheckSegments2(void) {
#if DEBUG_CURL_SEGMENT
 
#define DEBUG_TMP_N_CURL 50
 
  unsigned nSeg = m_segmentList.length();
  unsigned nWire[DEBUG_TMP_N_CURL] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                       0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  for(unsigned i=0;i<nSeg;++i){
    if(m_segmentList[i]->list().length() < DEBUG_TMP_N_CURL)
      ++(nWire[m_segmentList[i]->list().length()]);
    else
      ++(nWire[DEBUG_TMP_N_CURL-1]);
  }
  std::ifstream fin("tmp.wire.data");
  unsigned nTotalWire[DEBUG_TMP_N_CURL] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
                        0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
  if(fin){
    for(int i=0;i<DEBUG_TMP_N_CURL;++i)fin >> nTotalWire[i];
  }else{
    std::cout << "(TCurlFinder) tmp.wire.data does not exist!" << std::endl;
  }
  for(int i=0;i<DEBUG_TMP_N_CURL;++i)nTotalWire[i] += nWire[i];
  std::ofstream fout("tmp.wire.data");
  if(fout){
    fout << nTotalWire[0];
    for(int i=1;i<DEBUG_TMP_N_CURL;++i)fout << " " << nTotalWire[i];
  }else{
    std::cout << "(TCurlFinder) tmp.wire.data can not be made!" << std::endl;
  }
#endif
  return;
}
 
#undef DEBUG_CURL_DUMP
#undef DEBUG_CURL_SEGMENT
#undef DEBUG_CURL_GNUPLOT
#undef DEBUG_CURL_MC