/home/bes3soft/bes3soft/Boss/7.0.2/dist/7.0.2/Reconstruction/TrkReco/TrkReco-00-08-59-patch4-slc6tag/src/TBuilderCurl.cxx

Go to the documentation of this file.

//-----------------------------------------------------------------------------
// $Id: TBuilderCurl.cxx,v 1.21 2012/05/28 05:16:29 maoh Exp $
//-----------------------------------------------------------------------------
// Filename : TBuilderCurl.cc
// Section  : Tracking
// Owner    : Yoshi Iwasaki
// Email    : yoshihito.iwasaki@kek.jp
//-----------------------------------------------------------------------------
// Description : A class to build a curl track.
//               See http://bsunsrv1.kek.jp/~yiwasaki/tracking/
//-----------------------------------------------------------------------------

//#ifdef TRKRECO_DEBUG_DETAIL
//#define TRKRECO_DEBUG
//#endif
#include "TrkReco/TBuilderCurl.h"
#include "TrkReco/TMLink.h"
#include "TrkReco/TTrack.h"
#include "TrackUtil/Lpav.h"
//#include "TrkReco/Lpav.h"
//#include "TrkReco/TSvdFinder.h"
//#include "TrkReco/TSvdHit.h"

#include "TrkReco/TMLine.h"
#include "TrkReco/TRobustLineFitter.h"

// Following 3 parameters : (0,0,0) is best!
// Other cases are for the development.
#define DEBUG_CURL_DUMP 0
#define DEBUG_CURL_GNUPLOT 0
#define DEBUG_CURL_MC 0

#if (DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
#include "TrkReco/TMDCWireHitMC.h"
#include "TrkReco/TTrackHEP.h"
//#include MDC_H
#include "MdcTables/MdcTables.h"
#include <algo.h>
#endif

#include "GaudiKernel/StatusCode.h"
#include "GaudiKernel/IInterface.h"
#include "GaudiKernel/Kernel.h"
#include "GaudiKernel/Service.h"
#include "GaudiKernel/ISvcLocator.h"
#include "GaudiKernel/SvcFactory.h"
#include "GaudiKernel/IDataProviderSvc.h"
#include "GaudiKernel/Bootstrap.h"
#include "GaudiKernel/MsgStream.h"
#include "GaudiKernel/SmartDataPtr.h"
#include "GaudiKernel/IMessageSvc.h"

#if (DEBUG_CURL_DUMP+DEBUG_CURL_GNUPLOT+DEBUG_CURL_MC)
//#include "tables/belletdf.h"
//#include "tables/bestdf.h"
static int debugMcFlag = 1;
#endif

TBuilderCurl::TBuilderCurl(const std::string & name)
: TBuilder0(name),
  _fitter("TBuilderCurl Fitter"){
 
  //jialk
  StatusCode scmgn = Gaudi::svcLocator()->service("MagneticFieldSvc",m_pmgnIMF); 
  if(scmgn!=StatusCode::SUCCESS) { 
    std::cout<< __FILE__<<"Unable to open Magnetic field service"<<std::endl;
  }
  
#if 0
//  if(m_param.ON_CORRECTION){
//    _fitter.sag(true);
//    _fitter.propagation(true);
//    _fitter.tof(true);
//  }
  if (m_param.ON_CORRECTION & 1) _fitter.sag(true);
  if (m_param.ON_CORRECTION & 2) _fitter.propagation(true);
  if (m_param.ON_CORRECTION & 4) _fitter.tof(true);
#endif
}

TBuilderCurl::~TBuilderCurl() {
  //  if(m_param.SVD_RECONSTRUCTION){
  //    delete m_svdFinder;
  //    delete m_svdAssociator;
  //  }
}

void 
TBuilderCurl::setParam(const TCurlFinderParameter &p) {
  m_param.Z_CUT = p.Z_CUT;
  m_param.Z_DIFF_FOR_LAST_ATTEND = p.Z_DIFF_FOR_LAST_ATTEND;
  m_param.SELECTOR_MAX_SIGMA = p.SELECTOR_MAX_SIGMA;
  m_param.SELECTOR_MAX_IMPACT = p.SELECTOR_MAX_IMPACT;
  m_param.SELECTOR_STRANGE_PZ = p.SELECTOR_STRANGE_PZ;
  m_param.SELECTOR_REPLACE_DZ = p.SELECTOR_REPLACE_DZ;
  m_param.SVD_RECONSTRUCTION = p.SVD_RECONSTRUCTION;
  m_param.MIN_SVD_ELECTRONS = p.MIN_SVD_ELECTRONS;
  m_param.ON_CORRECTION = p.ON_CORRECTION;
  m_param.CURL_VERSION = p.CURL_VERSION;
#if 1
//  if(m_param.ON_CORRECTION){
//    _fitter.sag(true);
//    _fitter.propagation(true);
//    _fitter.tof(true);
//  }
  if (m_param.ON_CORRECTION & 1) _fitter.sag(true);
  if (m_param.ON_CORRECTION & 2) _fitter.propagation(true);
  if (m_param.ON_CORRECTION & 4) _fitter.tof(true);
#endif
  //  if(m_param.SVD_RECONSTRUCTION){
  //m_svdAssociator = new TSvdAssociator(-1.*(m_param.MIN_SVD_ELECTRONS),
  //                                     m_param.MIN_SVD_ELECTRONS);
  //}
}

// Stereo Finder For Curl Tracks : by jtanaka == from here ==

void 
TBuilderCurl::resetHelixFit(THelixFitter *fit) const {
  fit->fit2D(false);
  fit->sag(false);
  fit->propagation(false);
  fit->tof(false);
  fit->freeT0(false);
}

//.............................................
//...............Main Part.....................
//.............................................

bool
TBuilderCurl::buildStereo(TTrack & track,
                          double &dZ,
                          double &tanL) const {
/*  
 if(!(m_param.SVD_RECONSTRUCTION))return false;

#if 0
  double q = 1.;
  if(track.helix().kappa()<0.)q = -1.;
  double h[5], chisq;
  if(m_svdFinder->recTrk(track.helix().center().x(),
                         track.helix().center().y(),
                         fabs(track.helix().radius()),
                         q,
                         0.2,
                         h[0], h[1], h[2], h[3], h[4], chisq)){
#if 0
    std::cout << "SVD = "
         << h[0] << " "
         << h[1] << " "
         << h[2] << " "
         << h[3] << " "
         << h[4];
    std::cout << ", chi2 = " << chisq << " pt=" << 1./h[2] << std::endl;
#endif
    dZ   = h[3];
    tanL = h[4];
    return true;
  }else{
    //std::cout << "SVD Fail....." << std::endl;
    return false;
  }
#else
  Helix th(track.helix());
  th.pivot(HepPoint3D(0.,0.,0.));
  AList<TSvdHit> cand;
  if(m_svdAssociator->recTrk(th,dZ,tanL,0.5,-1.0,cand,0.5)){
    //std::cout << "SVD " << dZ << " " << tanL << std::endl;
    return true;
  }else{
    //std::cout << "SVD..." << std::endl;
    return false;
  }
#endif
 */
}

TTrack *
TBuilderCurl::buildStereo(TTrack & track, 
                          const AList<TMLink> & stereoList) const {
  return NULL;
}

TTrack *
TBuilderCurl::buildStereo(TTrack & track, 
                          const AList<TMLink> & stereoList,
                          const AList<TMLink> & allStereoList) const {
#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
#endif

  AList<TMLink> list = stereoList;
  unsigned nList = list.length();

  //...gets stereo wires from track.links
  for(unsigned i = 0, size = track.links().length(); i < size; ++i){
    unsigned superID = (track.links())[i]->wire()->superLayerId();
    if(superID == 0 || superID == 1 || superID == 5 || 
       superID == 6 || superID == 7 || superID == 8  ){
      int ok = 1;
      for (unsigned j = 0; j < nList; ++j) {
        TMLink * l = list[j];
        if(l->hit()->wire()->id() == (track.links())[i]->hit()->wire()->id()){
          ok = 0; break;
        }
      }
      if(ok == 1)list.append(((track.links())[i]));
    }
    // set LR 2
    (track.links())[i]->leftRight(2);  // returns left-right. 0:left, 1:right, 2:wire
  }
  for(unsigned i = 0, size = list.length(); i < size; ++i){    
    track._links.remove(list[i]);
    // set LR 2
    list[i]->leftRight(2);
  }

  if(list.length() < 2 || 
     list.length()+track.nLinks() < 5)return NULL;  
#if DEBUG_CURL_DUMP
  unsigned debug_stereo_counter1 = 0;
  for(unsigned i=0;i<track.links().length();++i){
    unsigned superID = (track.links())[i]->hit()->wire()->superLayerId();
    if(superID == 0 || superID == 1 ||
       superID == 5 || superID == 6 ||
       superID == 7 || superID == 8 )++debug_stereo_counter1;
  }
  std::cout << "(TBuilderCurl)Fitted Track:";
  std::cout << ", A# = " << track.links().length()-debug_stereo_counter1;
  std::cout << ", S# = " << debug_stereo_counter1 << "(==0)";
  std::cout << ", A+S# = " << track.links().length();
  std::cout << ", Cand Stereo Wires # = " << list.length() << std::endl;
  double debugChg = -1.;
  if(track.charge() > 0.)debugChg = 1.;
  if(debugChg > 0.)std::cout << "(TBuilderCurl)  Positive Track" << std::endl;
  else std::cout << "(TBuilderCurl)  Negative Track" << std::endl;
#endif

  //...calculates a circle.
  double xc2D;
  double yc2D;
  double r2D;
  AList<TMLink> tmpAxialList = track.links();
  bool err2D = fitWDD(xc2D,yc2D,r2D,tmpAxialList);

  //...using SVD   //....Liuqg
  //double dZSVD, tanLSVD;
  //bool initWithSVD = buildStereo(track, dZSVD, tanLSVD);

  //...sets arc and z pairs of each stereo hit wire.
  setArcZ(track,list);
  AList<TMLink> removeList;
  for(unsigned i = 0, size = list.length(); i < size; ++i){
    if(list[i]->arcZ(0).x() == -999. && list[i]->arcZ(1).x() == -999.)removeList.append(list[i]);
  }
  list.remove(removeList);
  if(list.length() < 2 || 
     list.length()+track.nLinks() < 5){ // stereo >=2 && axial >= 3
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)  Fail...few wires which can be set Arc-Z # = " 
         << list.length() << std::endl;
#endif    
    return NULL;
  }
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)  Cand Stereo Wires which can be set Arc-Z # = " 
         << list.length() << std::endl;
    plotArcZ(list,0.,0.,0);
#endif 
  
#if 0
  //...separates
  AList<TMLink> layer[18];
  for(unsigned i = 0, size = list.length(); i < size; ++i){
    layer[list[i]->wire()->layer()->axialStereoLayerId()].append(*list[i]);
  }

#if DEBUG_CURL_DUMP
  for(int i=0;i<18;++i){
    if(layer[i].length())
      std::cout << "(TBuilderCurl)  Cand Stereo Wires which can be set Arc-Z # = " 
           << layer[i].length()
           << " on " << i << " Layer." << std::endl;
  }
#endif
#endif

  // makes a line.
  AList<TMLink> goodL;
  AList<HepPoint3D> goodP;
  double minChi2 = 9999.;
  double goodA   = 9999.;
  double goodB   = 9999.;
  makeLine(track, list, allStereoList, goodL, 
           minChi2, goodA, goodB, goodP);
  HepAListDeleteAll(goodP);

#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    a = " << goodA << ", b = " << goodB 
            << " (after makeLine-function)" << std::endl;
#endif

  // refits
  static TRobustLineFitter rfitter("Can you work well?");
  TMLine newLine0(goodL);
  if(rfitter.fit(newLine0) != 0){
#if DEBUG_CURL_DUMP
      std::cout << "(TBuilderCurl)    Fail...linefitting...fail." << std::endl;
#endif
    return NULL;
  }

#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    a = " << newLine0.a() << ", b = " << newLine0.b() 
            << " (after robustline-fit)" << std::endl;
#endif

  //...appends at last chance
  unsigned nGoodL = goodL.length();
  list.remove(goodL);
  AList<TMLink> goodL2;
  if(m_param.CURL_VERSION == 0){ // default
    if(fabs(newLine0.b()) < 10.){
      appendPoints(list, goodL2,
                   newLine0.a(), newLine0.b(), track, m_param.Z_DIFF_FOR_LAST_ATTEND);
    }else{ // in case of bad result of robustLineFitter. (2001/04/04)
      appendPoints(list, goodL2,
                   goodA, goodB, track, m_param.Z_DIFF_FOR_LAST_ATTEND);    
    }
  }else{
    // same with b20010409_2122 at least
    appendPoints(list, goodL2,
                 newLine0.a(), newLine0.b(), track, m_param.Z_DIFF_FOR_LAST_ATTEND);
  }
  goodL.append(goodL2);
  TMLine newLine(goodL);

  //...refits
  if(rfitter.fit(newLine) != 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)    Fail...appending and re-fitting...fail." << std::endl;
#endif
    return NULL;
  }

#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    a = " << newLine.a() << ", b = " << newLine.b() 
            << " (after last-append + re-robustline-fit)" << std::endl;
#endif
  //...makes 3D tracks
  const AList<TMLink> &good = newLine.links();
  track.TTrackBase::append(good);
  if(!check(track)){
#if DEBUG_CURL_DUMP
      std::cout << "(TBuilderCurl)    Fail...checking wire numbers...fail." << std::endl;
#endif
    return NULL;
  }

  //...sets initial values
  Vector a = track.helix().a();
  if(err2D){
    double tmpA[3];
    double tmpQ = track.charge() > 0. ? 1. : -1.;
    tmpA[1] = fmod(atan2(tmpQ * yc2D,
                         tmpQ * xc2D)
                   + 4. * M_PI,
                   2. * M_PI);
   // tmpA[2] = tmpQ*Helix::ConstantAlpha/r2D;
   // tmpA[2] = tmpQ*333.564095/r2D;
    tmpA[2] = tmpQ*(333.564095/(-1000*(m_pmgnIMF->getReferField())))/r2D;
//cout<<"TBuilderCurl::tmpA[2]: "<<(333.564095/(-1000*(m_pmgnIMF->getReferField())))<<endl;
    tmpA[0] = xc2D/cos(tmpA[1]) - tmpQ*r2D;
    //std::cout << yc2D/sin(tmpA[1])-tmpQ*r2D << std::endl;
    //std::cout << a[0] << " -0- " << tmpA[0] << std::endl;
    //std::cout << a[1] << " -1- " << tmpA[1] << std::endl;
    //std::cout << a[2] << " -2- " << tmpA[2] << std::endl;
    a[0] = tmpA[0];
    a[1] = tmpA[1];
    a[2] = tmpA[2];
  }
  a[3] = newLine.b();
  a[4] = newLine.a();
//Liuqg
/*  if(initWithSVD){ // use SVD initial values if possible.
    a[3] = dZSVD;
    a[4] = tanLSVD;
  }
*/
  if(m_param.CURL_VERSION == 0){
    if(fabs(a[3]) > 10. && fabs(goodB) < 10.){ // 50cm --> 10cm @ 2001/04/04
      // use initial values when results of RobustFit is bad.
      a[3] = goodB;
      a[4] = goodA;
    }
  }else{
    if(fabs(a[3]) > 50. && fabs(goodB) < 50.){
      a[3] = goodB;
      a[4] = goodA;
      track.TTrackBase::remove(goodL2);
    }
  }

  track._helix->a(a);

#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    Created Line: y = " << newLine.a() 
       << " * x + " << newLine.b() 
       << ", size = " << good.length() << std::endl;
  plotArcZ(const_cast< AList<TMLink>& >(good), newLine.a(), newLine.b(), 0);
#endif

  /* std::cout << "1st : "
       << track.helix().a()[0] << " " << track.helix().a()[1] << " "
       << track.helix().a()[2] << " " << track.helix().a()[3] << " "
       << track.helix().a()[4] << std::endl; */

#if 1
//  if(m_param.ON_CORRECTION){
//    _fitter.sag();
//    _fitter.propagation();
//    _fitter.tof();
//  }
  if (m_param.ON_CORRECTION & 1) _fitter.sag(true);
  if (m_param.ON_CORRECTION & 2) _fitter.propagation(true);
  if (m_param.ON_CORRECTION & 4) _fitter.tof(true);
#endif

  //...fits
  AList<TMLink> bad;
  int err = _fitter.fit(track);
  if (err < 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail fitting(0)...error code = " << err << std::endl;
#endif
    return NULL;
  } else if ( fabs(track.helix().a()[3]) > fabs(a[3]) &&
              (fabs(track.helix().a()[3]) > 50.  || // 50cm
               fabs(track.helix().a()[2]) > 100. || // 0.01GeV
               fabs(track.helix().a()[2]) < 0.1) ){ // 10GeV
    // in strange case, set "correction" of wires OFF
    // and then fit with the initial values.
    if(m_param.ON_CORRECTION){
      if (m_param.ON_CORRECTION & 1) _fitter.sag(true);
      if (m_param.ON_CORRECTION & 2) _fitter.propagation(true);
      if (m_param.ON_CORRECTION & 4) _fitter.tof(true);
      //_fitter.sag();
      //_fitter.propagation();
      //_fitter.tof();
      track._helix->a(a);
      //std::cout << "ON --> OFF" << std::endl;
    }else{
      return NULL;
    }
  }

  /* std::cout << "2nd : "
       << track.helix().a()[0] << " " << track.helix().a()[1] << " "
       << track.helix().a()[2] << " " << track.helix().a()[3] << " "
       << track.helix().a()[4] << std::endl; */

  track.refine(bad, m_param.SELECTOR_MAX_SIGMA * 100.);
  if (!check(track)){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail checking(1)..." << std::endl;
#endif
    return NULL;
  }
  err = _fitter.fit(track);
  if (err < 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail fitting(1)...error code = " << err << std::endl;
#endif
    return NULL;
  }
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA * 10.);
  if (!check(track)){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail checking(2)..." << std::endl;
#endif
    return NULL;
  }
  err = _fitter.fit(track);
  if (err < 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail fitting(2)...error code = " << err << std::endl;
#endif
    return NULL;
  }
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA);
  if (!check(track)){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail checking(3)..." << std::endl;
#endif
    return NULL;
  }
  err = _fitter.fit(track);
  if (err < 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail fitting(3)...error code = " << err << std::endl;
#endif
    return NULL;
  }

#if 0
  for(int i=0;i<track.links().length();++i){
    if(track.links()[i]->pull() > 36.){
      std::cout << track.links()[i]->wire()->id() 
           << " :+: " << track.links()[i]->pull() << std::endl;
    }
    if(track.links()[i]->hit()->state() & WireHitInvalidForFit)
      std::cout << "Not Valid For Fit!" << std::endl;
    //if(!(track.links()[i]->hit()->state() & WireHitFittingValid))
    //std::cout << "No-Valid For Fit!" << std::endl;
  }
#endif
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA);
  if (!check(track)){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail checking(4)..." << std::endl;
#endif
    return NULL;
  }
  err = _fitter.fit(track);
  if (err < 0){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)     Fail fitting(4)...error code = " << err << std::endl;
#endif
    return NULL;
  }

#if 0    
  for(int i=0;i<track.links().length();++i){
    if(track.links()[i]->pull() > 36.){
      std::cout << track.links()[i]->wire()->id() 
           << " :*: " << track.links()[i]->pull() << std::endl;
    }
    if(track.links()[i]->hit()->state() & WireHitInvalidForFit)
      std::cout << "Not Valid For Fit!" << std::endl;
    //if(!(track.links()[i]->hit()->state() & WireHitFittingValid))
    //std::cout << "No-Valid For Fit!" << std::endl;
  }
#endif

  /* std::cout << "3rd : "
       << track.helix().a()[0] << " " << track.helix().a()[1] << " "
       << track.helix().a()[2] << " " << track.helix().a()[3] << " "
       << track.helix().a()[4] << std::endl; */

  //...tests
  if (track.nLinks() < 5){ // axial + stereo >= 5
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)      Success fitting, but pre-selection...fail." << std::endl;
#endif
    return NULL;
  }
  if(fabs(track.impact()) > m_param.SELECTOR_MAX_IMPACT ||
     track.pt() < 0.005 || // Pt >= 5MeV
     fabs(track.pz()) > m_param.SELECTOR_STRANGE_PZ){
#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)      Success fitting, but selection...fail." << std::endl;
    std::cout << "(TBuilderCurl)                           impact = " << track.impact() << std::endl;
    std::cout << "(TBuilderCurl)                           pt     = " << track.pt() << std::endl;
    std::cout << "(TBuilderCurl)                           pz     = " << track.pz() << std::endl;
    std::cout << "(TBuilderCurl)                           dz     = " << track.helix().a()[3] << std::endl;
#endif  
    return NULL;
  }

  //...replaces init helix if dz is bad.
  if(fabs(track.helix().a()[3]) > m_param.SELECTOR_REPLACE_DZ &&
     fabs(a[3]) < fabs(track.helix().a()[3])){
    track._helix->a(a);
  }
#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)      Success Build Stereo!!!" << std::endl;
#endif
  return & track;
}

//.............................................
//.............Set Arc and Z Part..............
//.............................................

void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &list) const {
  if(track.nLinks() < 4){
    track.stereoHitForCurl(list);
    return;
  }
//Liuqg  
  AList<TMLink> alayer[5];
  AList<TMLink> slayer[6];
  for(unsigned i=0,size=track.nLinks();i<size;++i){
    unsigned id = (track.links())[i]->wire()->superLayerId();
    if(id == 2)alayer[0].append((track.links())[i]);
    else if(id ==  3)alayer[1].append((track.links())[i]);
    else if(id ==  4)alayer[2].append((track.links())[i]);
    else if(id ==  9)alayer[3].append((track.links())[i]);
    else if(id ==  10)alayer[4].append((track.links())[i]);
  }

  for(unsigned i=0,size=list.length();i<size;++i){
    unsigned id = list[i]->wire()->superLayerId();
    if(id == 0)slayer[0].append(list[i]);
    else if(id == 1)slayer[1].append(list[i]);
    else if(id == 5)slayer[2].append(list[i]);
    else if(id == 6)slayer[3].append(list[i]);
    else if(id == 7)slayer[4].append(list[i]);
    else if(id == 8)slayer[5].append(list[i]);
  }
#if 0
  std::cout << "Stereo = " 
       << slayer[0].length() << " "
       << slayer[1].length() << " "
       << slayer[2].length() << " "
       << slayer[3].length() << " "
       << slayer[4].length() << std::endl;
  std::cout << "Axial  = " 
       << alayer[0].length() << " "
       << alayer[1].length() << " "
       << alayer[2].length() << " "
       << alayer[3].length() << " "
       << alayer[4].length() << " "
       << alayer[5].length() << std::endl;
#endif
  unsigned ip = 0;
  if(slayer[0].length() >= 1){
    if(alayer[0].length()+alayer[1].length() >= 4){
      setArcZ(track,slayer[0],alayer[0],alayer[1],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length() >= 4){
      setArcZ(track,slayer[0],alayer[0],alayer[1],
              alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length() >= 4){
       setArcZ(track,slayer[0],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[0],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);
    }
  }
  if(slayer[1].length() >= 1){
    if(alayer[0].length()+alayer[1].length() >= 4){
      setArcZ(track,slayer[1],alayer[0],alayer[1],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length() >= 4){
      setArcZ(track,slayer[1],alayer[0],alayer[1],
              alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length() >= 4){
      setArcZ(track,slayer[1],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[1],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);
    }
  }
  if(slayer[2].length() >= 1){
    if(alayer[1].length()+alayer[2].length() >= 4){
      setArcZ(track,slayer[2],alayer[1],alayer[2],ip);
    }
    else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length() >= 4){
      setArcZ(track,slayer[2],alayer[0],alayer[1],alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length() >= 4){
      setArcZ(track,slayer[2],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[2],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);

    }
  }
  if(slayer[3].length() >= 1){
    if(alayer[1].length()+alayer[2].length() >= 4){
      setArcZ(track,slayer[3],alayer[1],alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length() >= 4){
      setArcZ(track,slayer[3],alayer[0],alayer[1],
              alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()
              >= 4){
      setArcZ(track,slayer[3],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[3],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);

    }
  }
  if(slayer[4].length() >= 1){
    if(alayer[1].length()+alayer[2].length()
              >= 4){
      setArcZ(track,slayer[4],alayer[1],alayer[2],
             ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length() >= 4){
      setArcZ(track,slayer[4],alayer[0],alayer[1],
              alayer[2],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()
              >= 4){
      setArcZ(track,slayer[4],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[4],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);
    }
  }
  if(slayer[5].length() >= 1){
    if(alayer[1].length()+alayer[2].length()
              >= 4){
      setArcZ(track,slayer[5],alayer[1],alayer[2],
             ip);
    }else if(alayer[1].length()+alayer[2].length()+
             alayer[3].length() >= 4){
      setArcZ(track,slayer[5],alayer[1],alayer[2],
              alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()
              >= 4){
      setArcZ(track,slayer[5],alayer[0],alayer[1],
              alayer[2],alayer[3],ip);
    }else if(alayer[0].length()+alayer[1].length()+
             alayer[2].length()+alayer[3].length()+
             alayer[4].length() >= 4){
      setArcZ(track,slayer[5],alayer[0],alayer[1],
              alayer[2],alayer[3],alayer[4],ip);
    }
  }
}
//.............................................
//.............Append Points(last).............
//.............................................

unsigned
TBuilderCurl::appendPoints(AList<TMLink> &list, AList<TMLink> &line,
                           double a, double b, TTrack &track, double z_cut) const {
  unsigned size = list.length();
  if(size == 0)return 0;
  unsigned counter(0);
  for(unsigned i=0;i<size;++i){
    for(unsigned j=0;j<4;++j){
      if(j <= 1 && list[i]->arcZ(j).x() == -999.)continue;
      else if(j > 1 && list[i]->arcZ(j).x() == -999.)break;
      double y = a*list[i]->arcZ(j).x()+b;
      if(fabs(y-list[i]->arcZ(j).y()) < z_cut){
        list[i]->position(list[i]->arcZ(j));
        line.append(list[i]);
        ++counter;
        break;
      }
    }
  }
  return counter;
}

//.............................................
//..............Line Fit Part..................
//.............................................

int
doLineFit(AList<TMLink> &points, 
          double &m_a, double &m_b, double &chi2, double &nhits,
          int ipC = 1)   //Liuqg, IP Constraint
{
    m_a = m_b = nhits = 0.;
    chi2 = 1.e+10;
    unsigned n = points.length();
    double sum = double(n);
    double sumX = 0., sumY = 0., sumX2 = 0., sumXY = 0., sumY2 = 0.;
    for (unsigned i = 0; i < n; i++) {
      TMLink & l = * points[i];
      
      double x = l.position().x();
      double y = l.position().y();
      sumX  += x;
      sumY  += y;
      sumX2 += x * x;
      sumXY += x * y;
      sumY2 += y * y;
    }
    if(ipC != 0)sum += 1.0;// IP Constraint
    nhits = sum;
    
    double m_det = sum * sumX2 - sumX * sumX;
    if(m_det == 0. && sum != 2.){
      return -1;
    }else if(m_det == 0. && sum == 2.){
      double x0 = points[0]->position().x();
      double y0 = points[0]->position().y();
      double x1 = points[1]->position().x();
      double y1 = points[1]->position().y();
      if(x0 == x1)return -1;
      m_a = (y0-y1)/(x0-x1);
      m_b = -m_a*x1 + y1;
      chi2 = 0.;
      return 0;
    }
    chi2 = 0.;
    m_a = (sumXY * sum - sumX * sumY) / m_det;
    m_b = (sumX2 * sumY - sumX * sumXY) / m_det;

    for(unsigned i = 0; i < n; i++) {
      TMLink & l = * points[i];
      
      double x = l.position().x();
      double y = l.position().y();
      double d =  y-(x*m_a+m_b);
      chi2 += d*d;
    }

    return 0;
}

void
TBuilderCurl::fitLine(AList<TMLink> &tmpLine, double &min_chi2, 
                      double &good_a, double &good_b, 
                      AList<TMLink> &goodLine, AList<HepPoint3D> &goodPosition,
                      int &overCounter) const {
#if 0
  // OLD
  unsigned size = tmpLine.length();
  TMLine line(tmpLine);
  if(!(line.fit())){
    double chi2 = line.chi2()/(static_cast<double>(size-2));
    if(chi2 < min_chi2 && fabs(line.b()) < m_param.Z_CUT){
      min_chi2 = chi2;
      good_a   = line.a();
      good_b   = line.b();
      goodLine = tmpLine;
      HepAListDeleteAll(goodPosition);
      for(unsigned i=0;i<size;++i)
        goodPosition.append(new HepPoint3D(const_cast<HepPoint3D&>(tmpLine[i]->position())));
    }
  }
  ++overCounter;
#else
  unsigned size = tmpLine.length();
  double ta,tb,tc,tn;
  if(!(doLineFit(tmpLine,ta,tb,tc,tn,1)) && tn > 2.){ // with IP
    double chi2 = tc/(tn-2.);
    if(chi2 < min_chi2 && fabs(tb) < m_param.Z_CUT){
      min_chi2 = chi2;
      good_a   = ta;
      good_b   = tb;
      goodLine = tmpLine;
      HepAListDeleteAll(goodPosition);
      for(unsigned i=0;i<size;++i)
        goodPosition.append(new HepPoint3D(const_cast<HepPoint3D&>(tmpLine[i]->position())));
    }
  }
  ++overCounter;
#endif
}

//.............................................
//.................Check Part..................
//.............................................

unsigned 
TBuilderCurl::check(const TTrack &track) const {
  unsigned nAhits(0), nShits(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())++nShits;
    else ++nAhits;
  }
  if(nAhits >= 3 && nShits >= 2)return 1;// hard coding
  return 0;
}

//.............................................
//.............Checker(Debuger)................
//.............................................

int
checkBorder(AList<TMLink> &layer0, 
            AList<TMLink> &layer1, 
            AList<TMLink> &layer2){
  AList<TMLink> list = layer0;
  list.append(layer1);
  list.append(layer2);
  int size = list.length();
  if(size <= 1)return 0;
  int layerId = list[0]->hit()->wire()->layerId();
  int maxLocalId = 79;
  if(layerId >= 15)maxLocalId = 127;
  if(layerId >= 23)maxLocalId = 159;
  if(layerId >= 32)maxLocalId = 207;
  if(layerId >= 41)maxLocalId = 255;
  int HId = (int)(0.8*(maxLocalId+1));
  int LId = (int)(0.2*(maxLocalId+1));
  int low  = 0;
  int high = 0;
  for(int i=0;i>size;++i){
    if(list[i]->hit()->wire()->localId() < LId)low = 1;
    else if(list[i]->hit()->wire()->localId() > HId)high = 1;
    if(low == 1 && high == 1)return 1;
  }
  return 0;
}

int
checkBorder(AList<TMLink> &layer0, 
            AList<TMLink> &layer1, 
            AList<TMLink> &layer2, 
            AList<TMLink> &layer3){
  AList<TMLink> list = layer0;
  list.append(layer1);
  list.append(layer2);
  list.append(layer3);
  int size = list.length();
  if(size <= 1)return 0;
  int layerId = list[0]->hit()->wire()->layerId();
  int maxLocalId = 79;
  if(layerId >= 15)maxLocalId = 127;
  if(layerId >= 23)maxLocalId = 159;
  if(layerId >= 32)maxLocalId = 207;
  if(layerId >= 41)maxLocalId = 255;
  int HId = (int)(0.8*(maxLocalId+1));
  int LId = (int)(0.2*(maxLocalId+1));
  int low  = 0;
  int high = 0;
  for(int i=0;i>size;++i){
    if(list[i]->hit()->wire()->localId() < LId)low = 1;
    else if(list[i]->hit()->wire()->localId() > HId)high = 1;
    if(low == 1 && high == 1)return 1;
  }
  return 0;
}

int
offsetBorder(TMLink *l){
  int layerId = l->hit()->wire()->layerId();
  int maxLocalId = 79;
  if(layerId >= 15)maxLocalId = 127;
  if(layerId >= 23)maxLocalId = 159;
  if(layerId >= 32)maxLocalId = 207;
  if(layerId >= 41)maxLocalId = 255;
  return maxLocalId+1;
}

void
makeList(AList<TMLink> &layer, AList<TMLink> &list, double q, int border, int checkB, TMLink *layer0){
  int n = layer.length();
  if(checkB == 0){
    for(int i=0;i<n;++i){
      if(q < 0.){
        if(layer[i]->hit()->wire()->localId() >= border)list.append(layer[i]);
      }else{
        if(layer[i]->hit()->wire()->localId() <= border)list.append(layer[i]);
      }
    }
  }else{
    //difficult!! --> puls offset
    int offset = offsetBorder(layer0);
    if(border*2 < offset)border += offset;
    for(int i=0;i<n;++i){
      int lId = layer[i]->hit()->wire()->localId();
      if(lId*2 < offset)lId += offset;
      if(q < 0.){
        if(lId >= border)list.append(layer[i]);
      }else{
        if(lId <= border)list.append(layer[i]);
      }
    }
  }
}

int
TBuilderCurl::sortByLocalId(AList<TMLink> &list) const{
  int size = list.length();
  if(size <= 1)return 0;
  int layerId = list[0]->hit()->wire()->layerId();
  int maxLocalId;
/*  if(layerId < 15)maxLocalId = 79;
  else if(layerId >= 15)maxLocalId = 127;
  else if(layerId >= 23)maxLocalId = 159;
  else if(layerId >= 32)maxLocalId = 207;
  else if(layerId >= 41)maxLocalId = 255;
*/
//Liuqg in this class, the parameter of superlayer changed to layerid.
  if(layerId == 0)maxLocalId = 39;
  else if(layerId == 1)maxLocalId = 43;
  else if(layerId == 2)maxLocalId = 47;
  else if(layerId == 3)maxLocalId = 55;
  else if(layerId == 4)maxLocalId = 63;
  else if(layerId == 5)maxLocalId = 71;
  else if(layerId < 8)maxLocalId = 79;
  else if(layerId < 24)maxLocalId = 159;
  else if(layerId < 28)maxLocalId = 175;
  else if(layerId < 32)maxLocalId = 207;
  else if(layerId < 43)maxLocalId = 239;
  int flag = 0;
  for(int i=0;i<size;++i){
    if(list[i]->hit()->wire()->localId() == 0 ||
       list[i]->hit()->wire()->localId() == 1 ||
       list[i]->hit()->wire()->localId() == maxLocalId-1 ||
       list[i]->hit()->wire()->localId() == maxLocalId){
      flag = 1;
      break;
    }
  }
  if(flag == 0)return 0;
  int maxDif = (int)(0.5*(maxLocalId+1));
  AList<TMLink> fList; //former
  AList<TMLink> lList; //later
  for(int i=0;i<size;++i){
    if(list[i]->hit()->wire()->localId() < maxDif)
      lList.append(list[i]);
    else
      fList.append(list[i]);
  }
  list.removeAll();
  list.append(fList);
  list.append(lList);
  if(fList.length() >= 1 &&
     lList.length() >= 1)return 1;
  else return 0;
}

//..................................................
//..................................................
//....................MC............................
//..................................................
//..................................................

TTrack *
TBuilderCurl::buildStereoMC(TTrack & track, const AList<TMLink> & stereoList) const {
#if DEBUG_CURL_MC
  AList<TMLink> list = stereoList;

  HepPoint3D center = track.helix().center();
  double r  = fabs(track.helix().curv());
  for(unsigned i = 0, size = list.length();
      i < size; ++i){
    HepPoint3D point((list[i]->hit()->mc()->datcdc()->m_xin+
                   list[i]->hit()->mc()->datcdc()->m_xout)*0.5,
                  (list[i]->hit()->mc()->datcdc()->m_yin+
                   list[i]->hit()->mc()->datcdc()->m_yout)*0.5,
                  0.);
    double cosdPhi = - center.dot((point - center).unit()) / center.mag();
    double dPhi;
    if(fabs(cosdPhi) < 1.0){
      dPhi = acos(cosdPhi);
    }else if(cosdPhi >= 1.0){
      dPhi = 0.;
    }else{
      dPhi = M_PI;
    }
    list[i]->position(HepPoint3D(r*dPhi,
                              (list[i]->hit()->mc()->datcdc()->m_zin+
                               list[i]->hit()->mc()->datcdc()->m_zout)*0.5,
                              0.));
    /* std::cout << list[i]->wire()->id() << ": " 
         << point.x() << " "
         << point.y() << " "
         << (list[i]->hit()->mc()->datcdc()->m_zin+
         list[i]->hit()->mc()->datcdc()->m_zout)*0.5 << std::endl; */
  }
  //std::cout << "A# = " << track.links().length() << ", S# = " << list.length() << std::endl;
  double xc2D;
  double yc2D;
  double r2D;
  bool err2D = fitWDD(xc2D,yc2D,r2D,track.links()); // axial only

  TLine0 newLine(list);
  if(newLine.fit() != 0) return NULL;
  const AList<TMLink> &good = newLine.links();
  track.append(good);
  Vector a = track.helix().a();
  a[3] = newLine.b();
  a[4] = newLine.a();
  if(err2D){
    double tmpA[3];
    double tmpQ = track.charge() > 0. ? 1. : -1.;
    tmpA[1] = fmod(atan2(tmpQ * yc2D,
                         tmpQ * xc2D)
                   + 4. * M_PI,
                   2. * M_PI);
   // tmpA[2] = tmpQ*Helix::ConstantAlpha/r2D;
   // tmpA[2] = tmpQ*333.564095/r2D;
    tmpA[2] =tmpQ*(333.564095/(-1000*(m_pmgnIMF->getReferField())))/r2D;
    tmpA[0] = xc2D/cos(tmpA[1]) - tmpQ*r2D;
    a[0] = tmpA[0];
    a[1] = tmpA[1];
    a[2] = tmpA[2];
  }
  track._helix->a(a);
#if 0
  std::cout << track.helix().a()[0] << " " << track.helix().a()[1] << " "
       << track.helix().a()[2] << " " << track.helix().a()[3] << " "
       << track.helix().a()[4] << std::endl;
#endif
  //...fits
  AList<TMLink> bad;
  int err = _fitter.fit(track);
  if (err < m_param.ERROR_FOR_HELIX_FIT)return NULL;
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA * 100.);
  err = _fitter.fit(track);
  if (err < m_param.ERROR_FOR_HELIX_FIT)return NULL;
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA * 10.);
  err = _fitter.fit(track);
  if (err < m_param.ERROR_FOR_HELIX_FIT)return NULL;
  track.refine(bad, m_param.SELECTOR_MAX_SIGMA);
  err = _fitter.fit(track);
  if (err < m_param.ERROR_FOR_HELIX_FIT)return NULL;
#if 0
  std::cout << track.helix().a()[0] << " " << track.helix().a()[1] << " "
       << track.helix().a()[2] << " " << track.helix().a()[3] << " "
       << track.helix().a()[4] << std::endl;
#endif
  //track._helix->a(a); // re-write

  //...checks
  if (track.nLinks() < m_param.SELECTOR_NLINKS)return NULL;
  if (fabs(track.impact()) > m_param.SELECTOR_MAX_IMPACT ||
      track.pt() < m_param.SELECTOR_MIN_PT)return NULL;
  return & track;
#else
  return NULL;
#endif
}

//..................................................
//....................Plot..........................
//..................................................

void
TBuilderCurl::plotArcZ(AList<TMLink> &tmpLine,
                       double a, double b, const int flag) const {
#if DEBUG_CURL_GNUPLOT
  //#if 1
  if(a == 9999. || b == 9999.){
    a = 0.;
    b = 0.;
  }
  int nCounter = 0;
  double gmaxX = 0. ,gminX = 0.;
  double gmaxY = 0. ,gminY = 0.;
  FILE *gnuplot, *data;
  if((data = fopen("you_can_remove_this.dat","w")) != NULL){
    for(int ii=0;ii<tmpLine.length();ii++){
      ++nCounter;
      fprintf(data,"%lf, %lf\n",
              tmpLine[ii]->position().x(),
              tmpLine[ii]->position().y());
      if(flag)std::cout << " Wire ID = " << tmpLine[ii]->hit()->wire()->id() 
                   << " Arc = " << tmpLine[ii]->position().x()
                   << " Z = " << tmpLine[ii]->position().y();
      //if(flag && !debugMcFlag)std::cout << std::endl;
      //if(flag && debugMcFlag){
      //std::cout << " Z(true) = " << (tmpLine[ii]->hit()->mc()->datcdc()->m_zin+
      //                          tmpLine[ii]->hit()->mc()->datcdc()->m_zout)*0.5;
      //std::cout << " HepTrackID = " << tmpLine[ii]->hit()->mc()->hep()->id() << std::endl;
      //}
      if(gmaxX < tmpLine[ii]->position().x())
        gmaxX = tmpLine[ii]->position().x();
      if(gminX > tmpLine[ii]->position().x())
        gminX = tmpLine[ii]->position().x();
      if(gmaxY < tmpLine[ii]->position().y())
        gmaxY = tmpLine[ii]->position().y();
      if(gminY > tmpLine[ii]->position().y())
        gminY = tmpLine[ii]->position().y();
    }
    fclose(data);
  }
  if((data = fopen("you_can_remove_this.dat2","w")) != NULL){
#if 0
    if(debugMcFlag){
      for(int ii=0;ii<tmpLine.length();ii++){
        double z = (tmpLine[ii]->hit()->mc()->datcdc()->m_zin+
                    tmpLine[ii]->hit()->mc()->datcdc()->m_zout)*0.5;
        if(tmpLine[ii]->arcZ(0).x() != -999.){
          ++nCounter;
          fprintf(data,"%lf, %lf\n",tmpLine[ii]->arcZ(0).x(),z);
          if(gmaxX < tmpLine[ii]->arcZ(0).x())
            gmaxX = tmpLine[ii]->arcZ(0).x();
          if(gminX > tmpLine[ii]->arcZ(0).x())
            gminX = tmpLine[ii]->arcZ(0).x();
          if(gmaxY < z)
            gmaxY = z;
          if(gminY > z)
            gminY = z;
        }
      }
    }
#endif
    fclose(data);
  }
  if((data = fopen("you_can_remove_this.dat3","w")) != NULL){
    for(int ii=0;ii<tmpLine.length();ii++){
      for(int jj=0;jj<4;++jj){
        if(tmpLine[ii]->arcZ(jj).x() != -999.){
          ++nCounter;
          fprintf(data,"%lf, %lf\n",
                  tmpLine[ii]->arcZ(jj).x(),
                  tmpLine[ii]->arcZ(jj).y());
          if(gmaxX < tmpLine[ii]->arcZ(jj).x())
            gmaxX = tmpLine[ii]->arcZ(jj).x();
          if(gminX > tmpLine[ii]->arcZ(jj).x())
            gminX = tmpLine[ii]->arcZ(jj).x();
          if(gmaxY < tmpLine[ii]->arcZ(jj).y())
            gmaxY = tmpLine[ii]->arcZ(jj).y();
          if(gminY > tmpLine[ii]->arcZ(jj).y())
            gminY = tmpLine[ii]->arcZ(jj).y();
        }else{
          break;
        }
      }
    }
    fclose(data);
  }
  if(gmaxX < 0.)gmaxX = 0.;if(gminX > 0.)gminX = 0.;
  if(gmaxY < 0.)gmaxY = 0.;if(gminY > 0.)gminY = 0.;
  if(nCounter && (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);
    if(a == 0. && b == 0.){
      fprintf(gnuplot,"plot \"you_can_remove_this.dat3\", \"you_can_remove_this.dat\", \"you_can_remove_this.dat2\" \n");
    }else{
      fprintf(gnuplot,"plot \"you_can_remove_this.dat3\", \"you_can_remove_this.dat\", \"you_can_remove_this.dat2\", %lf*x+%lf \n", a, b);
    }
    fflush(gnuplot);
    char tmp[8];
    gets(tmp);
    pclose(gnuplot);
  }
#endif
  return;
}

// 
//
//

unsigned
findMaxLocalId(unsigned layerId){
/*  unsigned maxLocalId = 79;
  if(superLayerId == 3)maxLocalId = 127;
  else if(superLayerId == 5)maxLocalId = 159;
  else if(superLayerId == 7)maxLocalId = 207;
  else if(superLayerId == 9)maxLocalId = 255;
  return maxLocalId;
*/
//changed to BESiii, Liuqg
  unsigned maxLocalId = 39;
  if(layerId == 1)maxLocalId = 43;
  else if(layerId == 2)maxLocalId = 47;
  else if(layerId == 3)maxLocalId = 55;
  else if(layerId == 4)maxLocalId = 63;
  else if(layerId == 5)maxLocalId = 71;
  else if(layerId == 6 || layerId == 7 )maxLocalId = 79;
  else if(layerId == 20 || layerId == 21 || layerId == 22 || layerId == 23)maxLocalId = 159;
  else if(layerId == 24 || layerId == 25 || layerId == 26 || layerId == 27)maxLocalId = 175;
  else if(layerId == 28 || layerId == 29 || layerId == 30 || layerId == 31)maxLocalId = 207;
  else if(layerId == 32 || layerId == 33 || layerId == 34 || layerId == 35)maxLocalId = 239;

  return maxLocalId;
}

unsigned
isIsolation(unsigned localId, 
              unsigned maxLocalId,
              unsigned layerId,
              int lr,
              const AList<TMLink> &allStereoList){
  unsigned findId;
  if(lr == 1){ // R : ox, Dose a wire exist at "x"?
    findId = maxLocalId;
    if(localId != 0)findId = localId-1;
  }else if(lr == -1){ // L : xo, Dose a wire exist at "x"?
    findId = 0;
    if(localId != maxLocalId)findId = localId+1;
  }else{
    return 1;
  }
  unsigned isolate = 1;
  for(int i=0;i<allStereoList.length();++i){
    if(allStereoList[i]->wire()->layerId() == layerId &&
       allStereoList[i]->wire()->localId() == findId){
      isolate = 0;
      break;     
    }
  }
  return isolate;
}

void
findTwoHits(AList<TMLink> &twoOnLayer,
            const AList<TMLink> &hitsOnLayer,
            const AList<TMLink> &allStereoList){
  //...finds "two" seq. and isolated hits.
  //...and then sets LR for selected hits.
  if(hitsOnLayer.length() == 0 ||
     hitsOnLayer.length() > 3)return;
  twoOnLayer.removeAll();
  if(hitsOnLayer.length() == 1){
    if(hitsOnLayer[0]->wire()->superLayerId() == 0)return; //origin is == 1, Liuqg 060921
    unsigned maxLocalId = findMaxLocalId(hitsOnLayer[0]->wire()->layerId());
    unsigned R = isIsolation(hitsOnLayer[0]->wire()->localId(),
                             maxLocalId,
                             hitsOnLayer[0]->wire()->layerId(),1,allStereoList);
    unsigned L = isIsolation(hitsOnLayer[0]->wire()->localId(),
                             maxLocalId,
                             hitsOnLayer[0]->wire()->layerId(),-1,allStereoList);
    if(R == 1 && L == 0){
      unsigned nextLocalId = hitsOnLayer[0]->wire()->localIdForPlus()+1;
      L = isIsolation(nextLocalId,
                      maxLocalId,
                      hitsOnLayer[0]->wire()->layerId(),-1,allStereoList);
      if(L == 1){ // xuox
        hitsOnLayer[0]->leftRight(1); // R
        hitsOnLayer[0]->position(hitsOnLayer[0]->arcZ(1));
        twoOnLayer.append(hitsOnLayer[0]);
      }
    }else if(R == 0 && L == 1){
      unsigned nextLocalId = hitsOnLayer[0]->wire()->localIdForMinus()+1;
      R = isIsolation(nextLocalId,
                      maxLocalId,
                      hitsOnLayer[0]->wire()->layerId(),1,allStereoList);
      if(R == 1){ // xoux
        hitsOnLayer[0]->leftRight(0); // L
        hitsOnLayer[0]->position(hitsOnLayer[0]->arcZ(0));
        twoOnLayer.append(hitsOnLayer[0]);
      }
    }   
  }
  if(hitsOnLayer.length() == 2){
    if(hitsOnLayer[0]->wire()->localIdForPlus()+1 ==
       hitsOnLayer[1]->wire()->localId()){
      unsigned maxLocalId = findMaxLocalId(hitsOnLayer[0]->wire()->layerId());
      unsigned R = isIsolation(hitsOnLayer[0]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[0]->wire()->layerId(),1,allStereoList);
      unsigned L = isIsolation(hitsOnLayer[1]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[1]->wire()->layerId(),-1,allStereoList);
      if(R == 1 && L == 1){ // xoox
        hitsOnLayer[0]->leftRight(1); // R
        hitsOnLayer[0]->position(hitsOnLayer[0]->arcZ(1));
        hitsOnLayer[1]->leftRight(0); // L
        hitsOnLayer[1]->position(hitsOnLayer[1]->arcZ(0));
        twoOnLayer.append(hitsOnLayer[0]);
        twoOnLayer.append(hitsOnLayer[1]);
      }
    }   
  }
  if(hitsOnLayer.length() == 3){
    if(hitsOnLayer[0]->wire()->localIdForPlus()+1 ==
       hitsOnLayer[1]->wire()->localId() &&
       hitsOnLayer[1]->wire()->localIdForPlus()+1 !=
       hitsOnLayer[2]->wire()->localId()){
      unsigned maxLocalId = findMaxLocalId(hitsOnLayer[0]->wire()->layerId());
      unsigned R = isIsolation(hitsOnLayer[0]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[0]->wire()->layerId(),1,allStereoList);
      unsigned L = isIsolation(hitsOnLayer[1]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[1]->wire()->layerId(),-1,allStereoList);
      if(R == 1 && L == 1){ // oxoox
        hitsOnLayer[0]->leftRight(1); // R
        hitsOnLayer[0]->position(hitsOnLayer[0]->arcZ(1));
        hitsOnLayer[1]->leftRight(0); // L
        hitsOnLayer[1]->position(hitsOnLayer[1]->arcZ(0));
        twoOnLayer.append(hitsOnLayer[0]);
        twoOnLayer.append(hitsOnLayer[1]);
      }
    }else if(hitsOnLayer[0]->wire()->localIdForPlus()+1 !=
             hitsOnLayer[1]->wire()->localId() &&
             hitsOnLayer[1]->wire()->localIdForPlus()+1 ==
             hitsOnLayer[2]->wire()->localId()){
      unsigned maxLocalId = findMaxLocalId(hitsOnLayer[1]->wire()->layerId());
      unsigned R = isIsolation(hitsOnLayer[1]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[1]->wire()->layerId(),1,allStereoList);
      unsigned L = isIsolation(hitsOnLayer[2]->wire()->localId(),
                               maxLocalId,
                               hitsOnLayer[2]->wire()->layerId(),-1,allStereoList);
      if(R == 1 && L == 1){ // xooxo
        hitsOnLayer[1]->leftRight(1); // R
        hitsOnLayer[1]->position(hitsOnLayer[1]->arcZ(1));
        hitsOnLayer[2]->leftRight(0); // L
        hitsOnLayer[2]->position(hitsOnLayer[2]->arcZ(0));
        twoOnLayer.append(hitsOnLayer[1]);
        twoOnLayer.append(hitsOnLayer[2]);
      }
    }
  }
  /* if(twoOnLayer.length() != 0){
    std::cout << "TWO " << twoOnLayer.length() << std::endl;
    } */
}

void
setLR(AList<TMLink> &hitsOnLayer, unsigned LR = 0){
  // LR = 0 : L
  //    = 1 : R
  for(unsigned i=0;i<hitsOnLayer.length();++i){
    if(LR == 0){
      hitsOnLayer[i]->leftRight(0); // L
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(0));
    }else{
      hitsOnLayer[i]->leftRight(1); // R
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(1));
    }
  }
}

bool
moveLR(AList<TMLink> &hitsOnLayer){
  unsigned nHits = hitsOnLayer.length();
  if(nHits == 0)return false;
  // ex) LLLL --> LLLR --> LLRR --> LRRR --> RRRR
  if(hitsOnLayer[nHits-1]->leftRight() == 1)return false; // All R
  for(unsigned i=0;i<nHits;++i){
    if(hitsOnLayer[i]->leftRight() == 0){ // L
      hitsOnLayer[i]->leftRight(1); // R
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(1));
      return true;
    }
  }
  return false;
}

void
selectGoodWires(const AList<TMLink> &allWires,
                AList<TMLink> &goodWires){
  goodWires.removeAll();
  for(int i=0;i<allWires.length();++i){
    if(allWires[i]->position().x() != -999.){
      goodWires.append(allWires[i]);
    }
  }
}

void
TBuilderCurl::fitLine2(const AList<TMLink> &tmpLine, double &min_chi2, 
                       double &good_a, double &good_b, 
                       AList<TMLink> &goodLine, AList<HepPoint3D> &goodPosition,
                       int &overCounter) const {
  AList<TMLink> goodWires;
  selectGoodWires(tmpLine,goodWires);
  if(goodWires.length() >= 3)
    fitLine(goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
}

void
calVirtualCircle(const TMLink &hit, const TTrack &track, const int LR,
                 HepPoint3D &center, double &radius){
  if(abs(LR) != 1)return;
  double Q = track.charge();
  int isOuter = 1;
  if(Q > 0. && LR == 1)isOuter = -1; // Inner
  else if(Q < 0. && LR == -1)isOuter = -1; // Inner
  radius = fabs(track.radius());
  center = track.helix().center();    
  HepPoint3D wire = hit.wire()->xyPosition();
  center.setZ(0.);
  wire.setZ(0.);
  // new center(virtual)
  center = wire + 
    (radius+((double)isOuter)*(hit.hit()->drift()))*(center-wire).unit();
}

void
moveLR(AList<TMLink> &hits,
       const AList<TMLink> &hitsOnLayerOrg,
       const TTrack &track){
  AList<TMLink> hitsOnLayer = hitsOnLayerOrg;
  hitsOnLayer.remove(hits);
  if(hitsOnLayer.length() == 0)return;

  unsigned nHits = hits.length();
  if(nHits == 0)return;

  //...finds "ref" from hits.
  // ex) LLLL|, LLL|R, LL|RR, L|RRR, |RRRR
  int LR = -1; // L
  TMLink ref;
  if(hits[nHits-1]->leftRight() == 1){ // All R
    LR = 1; // R
    ref = *hits[nHits-1];
  }
  for(unsigned i=0;i<nHits;++i){
    if(hits[i]->leftRight() == 0){ // L
      ref = *hits[i];
      break;
    }
  }

  HepPoint3D center;
  double radius;
  calVirtualCircle(ref,track,LR,center,radius);

  double Q = track.charge();
  for(int i=0;i<hitsOnLayer.length();++i){
    int isOuter = 1;
    if((hitsOnLayer[i]->wire()->xyPosition()-center).mag()-radius < 0.)isOuter = -1;
    if(Q > 0. && isOuter == 1){
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(0)); // L
      hitsOnLayer[i]->leftRight(0); // L
    }else if(Q > 0. && isOuter == -1){
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(1)); // R
      hitsOnLayer[i]->leftRight(1); // R
    }else if(Q < 0. && isOuter == 1){
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(1)); // R
      hitsOnLayer[i]->leftRight(1); // R
    }else if(Q < 0. && isOuter == -1){
      hitsOnLayer[i]->position(hitsOnLayer[i]->arcZ(0)); // L
      hitsOnLayer[i]->leftRight(0); // L
    }
  }
}

void
TBuilderCurl::makeLine(TTrack &track, 
                       AList<TMLink> &list, const AList<TMLink> &allStereoList, AList<TMLink> &goodLine, 
                       double &min_chi2, double &good_a, double &good_b,
                       AList<HepPoint3D> &goodPosition) const {
  if(list.length() == 0)return;

  AList<TMLink> layer[24]; //Liuqg, origin is 18.
  AList<TMLink> layerOrg[24]; //Liuqg, origin is 18.
  for(unsigned i = 0, size = list.length(); i < size; ++i){
    layer[list[i]->wire()->layer()->axialStereoLayerId()].append(*list[i]);
    layerOrg[list[i]->wire()->layer()->axialStereoLayerId()].append(*list[i]);
  }

  AList<TMLink> fixedWires[6]; // each superlayer    origin is 5,Liuqg 060920
  AList<TMLink> nonFixedWires[6]; // each superlayer    origin is 5,Liuqg 060920
  AList<TMLink> allFixedWires;
  for(unsigned i=0;i<24;++i){
    if(layer[i].length()){
      layer[i].sort(SortByWireId);
      sortByLocalId(layer[i]); // chiisai-jun but kyoukai fukin ha sukoshi kufuu shite iru.
      AList<TMLink> tmp;
      findTwoHits(tmp,layer[i],allStereoList);
      if(tmp.length()){
        layer[i].removeAll();
        allFixedWires.append(tmp);
//the following parameters have been changed, liuqg 060919
        if(i<4)fixedWires[0].append(tmp);
        else if(i<8)fixedWires[1].append(tmp);
        else if(i<12)fixedWires[2].append(tmp);
        else if(i<16)fixedWires[3].append(tmp);
        else if(i<20)fixedWires[4].append(tmp);
        else fixedWires[5].append(tmp);
      }else{
        if(i<4)nonFixedWires[0].append(layer[i]);
        else if(i<8)nonFixedWires[1].append(layer[i]);
        else if(i<12)nonFixedWires[2].append(layer[i]);
        else if(i<16)nonFixedWires[3].append(layer[i]);
        else if(i<20)nonFixedWires[4].append(layer[i]);
        else nonFixedWires[5].append(layer[i]);
      }      
    }
  }

#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    1st fixed & non-fixed wires selection:" << std::endl;
  std::cout << "(TBuilderCurl)    all fixed wires # = " << allFixedWires.length() << std::endl; 
  std::cout << "(TBuilderCurl)    fixed wires # = (" 
       << fixedWires[0].length() << ", "
       << fixedWires[1].length() << ", "
       << fixedWires[2].length() << ", "
       << fixedWires[3].length() << ", "
       << fixedWires[4].length() << ", "
       << fixedWires[5].length() << ")" << std::endl;
  std::cout << "(TBuilderCurl)    non fixed wires # = (" 
       << nonFixedWires[0].length() << ", "
       << nonFixedWires[1].length() << ", "
       << nonFixedWires[2].length() << ", "
       << nonFixedWires[3].length() << ", "
       << nonFixedWires[4].length() << ", "
       << nonFixedWires[5].length() << ")" << std::endl;
#if 0 /* in detail */ 
  for(unsigned i=0;i<allFixedWires.length();++i)
    std::cout << i << ": LocalID/LayerID = " << allFixedWires[i]->wire()->localId() 
         << "/" << allFixedWires[i]->wire()->layerId()
         << ", LR = " << allFixedWires[i]->leftRight() << std::endl;
#endif /* in detail */
#endif

  int createdLine = 0;
  int overCounter = 0;
  AList<TMLink> goodWires;
#if 1 /* fastest finder */
  if(allFixedWires.length() >= 5){
    selectGoodWires(allFixedWires,goodWires);
    if(goodWires.length() >= 5){
      fitLine(goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
      if(fabs(good_b) < 10.)createdLine = 1;
    }
  }
#endif /* fastest finder */
#if 1 /* faster finder */
  if(createdLine == 0){
    // Q > 0 : Outer = L, Inner = R
    // Q < 0 : Outer = R, Inner = L
    double Q = track.charge();
    unsigned isIncreased = 0;
    for(int sl=0;sl<6;++sl){        //origin is 5, Liuqg 060919
      if(fixedWires[sl].length()    >= 1 &&
         nonFixedWires[sl].length() >= 1){
        isIncreased = 1;
        unsigned bestNCorrectLR = 0;
        double bestR;
        HepPoint3D bestC;
        for(int i=0;i<fixedWires[sl].length();++i){
          int LR = fixedWires[sl][i]->leftRight() == 0 ? -1 : 1;
          HepPoint3D center;
          double radius;
          calVirtualCircle(*fixedWires[sl][i],track,LR,center,radius);
          unsigned nCorrectLR = 0;
          for(int j=0;j<fixedWires[sl].length();++j){
            if(i != j){
              int tmpIsOuter = 1;
              if((fixedWires[sl][j]->wire()->xyPosition()-center).mag()-radius < 0.)tmpIsOuter = -1;
              if(Q > 0. && tmpIsOuter == 1 && fixedWires[sl][j]->leftRight() == 0)++nCorrectLR;
              else if(Q > 0. && tmpIsOuter == -1 && fixedWires[sl][j]->leftRight() == 1)++nCorrectLR;
              else if(Q < 0. && tmpIsOuter ==  1 && fixedWires[sl][j]->leftRight() == 1)++nCorrectLR;
              else if(Q < 0. && tmpIsOuter == -1 && fixedWires[sl][j]->leftRight() == 0)++nCorrectLR;
            }
          }
          if(i == 0 || nCorrectLR > bestNCorrectLR){
            bestNCorrectLR = nCorrectLR;
            bestR = radius;
            bestC = center;
          }
          if(bestNCorrectLR == fixedWires[sl].length()-1)break;
        }
        for(int i=0;i<nonFixedWires[sl].length();++i){
          int isOuter = 1;
          if((nonFixedWires[sl][i]->wire()->xyPosition()-bestC).mag()-bestR < 0.)isOuter = -1;
          if(Q > 0. && isOuter == 1){
            nonFixedWires[sl][i]->position(nonFixedWires[sl][i]->arcZ(0)); // L
            nonFixedWires[sl][i]->leftRight(0); // L
          }else if(Q > 0. && isOuter == -1){
            nonFixedWires[sl][i]->position(nonFixedWires[sl][i]->arcZ(1)); // R
            nonFixedWires[sl][i]->leftRight(1); // R
          }else if(Q < 0. && isOuter == 1){
            nonFixedWires[sl][i]->position(nonFixedWires[sl][i]->arcZ(1)); // R
            nonFixedWires[sl][i]->leftRight(1); // R
          }else if(Q < 0. && isOuter == -1){
            nonFixedWires[sl][i]->position(nonFixedWires[sl][i]->arcZ(0)); // L
            nonFixedWires[sl][i]->leftRight(0); // L
          }
        }
        allFixedWires.append(nonFixedWires[sl]);
        fixedWires[sl].append(nonFixedWires[sl]);
        nonFixedWires[sl].removeAll();
      }
    }

#if DEBUG_CURL_DUMP
    std::cout << "(TBuilderCurl)    2nd fixed & non-fixed wires selection:" << std::endl;
    std::cout << "(TBuilderCurl)    all fixed wires # = " << allFixedWires.length() << std::endl; 
    std::cout << "(TBuilderCurl)    fixed wires # = (" 
         << fixedWires[0].length() << ", "
         << fixedWires[1].length() << ", "
         << fixedWires[2].length() << ", "
         << fixedWires[3].length() << ", "
         << fixedWires[4].length() << ", "
         << fixedWires[5].length() << ")" << std::endl;
    std::cout << "(TBuilderCurl)    non fixed wires # = (" 
         << nonFixedWires[0].length() << ", "
         << nonFixedWires[1].length() << ", "
         << nonFixedWires[2].length() << ", "
         << nonFixedWires[3].length() << ", "
         << nonFixedWires[4].length() << ", "
         << nonFixedWires[5].length() << ")" << std::endl;
#if 0 /* in detail */ 
    for(unsigned i=0;i<allFixedWires.length();++i)
      std::cout << i << ": LocalID/LayerID = " << allFixedWires[i]->wire()->localId() 
           << "/" << allFixedWires[i]->wire()->layerId()
           << ", LR = " << allFixedWires[i]->leftRight() << std::endl;
#endif /* in detail */
#endif

    if(isIncreased == 1 && allFixedWires.length() >= 5){
      selectGoodWires(allFixedWires,goodWires);
      if(goodWires.length() >= 5){
        fitLine(goodWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
        if(fabs(good_b) < 10.)createdLine = 1;
      }
    }
  }
#endif /* faster finder */
#if 1 /* slow finder */
  // 2000/1/27...very slow but unlike an infinity loop
  if(createdLine == 0){
    // nonFixed Wires
    int maxNonFixedLayerIndex[6] = { -1, -1, -1, -1, -1, -1 }; //origin is 5, Liuqg 060919
    int maxLength[6] = { 0, 0, 0, 0, 0, 0 }; //origin is 5, Liuqg 060919
    for(int l=0;l<24;++l){
      unsigned sl = 5; // superlayer id, changed by Liuqg
      if(l<4)sl = 0;
      else if(l<8)sl = 1;
      else if(l<12)sl = 2;
      else if(l<16)sl = 3;
      else if(l<20)sl = 4;
      layer[l].remove(fixedWires[sl]);
      setLR(layer[l]); // set All L
      if(layer[l].length() && layer[l].length() > maxLength[sl]){
        maxLength[sl] = layer[l].length();
        maxNonFixedLayerIndex[sl] = l;
      }
    }

    unsigned index = 0;
    unsigned nonFixedSuperLayerIndex[6] = { 0,0,0,0,0,0 };   //origin is 5, Liuqg 060919
    unsigned isIncreased = 0;
    for(unsigned i=0;i<6;++i){  //origin is 5, Liuqg 060919
      allFixedWires.append(nonFixedWires[i]);
      if(nonFixedWires[i].length()){
        isIncreased = 1;
        nonFixedSuperLayerIndex[index] = i;
        ++index;
      }
    }

    if(isIncreased){
      do{
        moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[0]]],
               nonFixedWires[nonFixedSuperLayerIndex[0]],
               track);
        if(index > 1){
          setLR(nonFixedWires[nonFixedSuperLayerIndex[1]]);
          do{
            moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[1]]],
                   nonFixedWires[nonFixedSuperLayerIndex[1]],
                   track);
            if(index > 2){
              setLR(nonFixedWires[nonFixedSuperLayerIndex[2]]);
              do{
                moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[2]]],
                       nonFixedWires[nonFixedSuperLayerIndex[2]],
                       track);
                if(index > 3){
                  setLR(nonFixedWires[nonFixedSuperLayerIndex[3]]);
                  do{
                    moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[3]]],
                           nonFixedWires[nonFixedSuperLayerIndex[3]],
                           track);
                    if(index > 4){
                      setLR(nonFixedWires[nonFixedSuperLayerIndex[4]]);
                      do{
                        moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[4]]],
                               nonFixedWires[nonFixedSuperLayerIndex[4]],
                               track);
                    if(index > 5){
                      setLR(nonFixedWires[nonFixedSuperLayerIndex[5]]);
                      do{
                        moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[5]]],
                               nonFixedWires[nonFixedSuperLayerIndex[5]],
                               track);
                        fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
                        if(overCounter>10000)goto kokohe;
                      }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[5]]]));
                    }else{
                      fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
                        if(overCounter>10000)goto kokohe;
                        }
                      }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[4]]]));
                    }else{
                      fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
                      if(overCounter>10000)goto kokohe;
                    }
                  }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[3]]]));
                }else{
                  fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
                  if(overCounter>10000)goto kokohe;
                }
              }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[2]]]));
            }else{
              fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
              if(overCounter>10000)goto kokohe;
            }
          }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[1]]]));
        }else{
          fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
          if(overCounter>10000)goto kokohe;
        }
      }while(moveLR(layer[maxNonFixedLayerIndex[nonFixedSuperLayerIndex[0]]]));
    kokohe:;
    }else if(allFixedWires.length() >= 3){
      fitLine2(allFixedWires, min_chi2, good_a, good_b, goodLine, goodPosition, overCounter);
    }
  }
#endif /* slow finder */
  for(unsigned i = 0, size = goodLine.length(); i < size; ++i){
    goodLine[i]->position(*(goodPosition[i]));
  }
#if DEBUG_CURL_DUMP
  std::cout << "(TBuilderCurl)    make a line from All SuperLayers." << std::endl;
  plotArcZ(goodLine, good_a, good_b, 0);
#endif
}

bool
TBuilderCurl::fitWDD(double &xc, double &yc, double &r,
                     AList<TMLink> &list) const {
  if(list.length() <= 3)return false;
  Lpav circle;
  // MDC
  for(int i=0;i<list.length();++i){
    circle.add_point(list[i]->wire()->xyPosition().x(),     
                     list[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;
  xc = circle.center()[0];
  yc = circle.center()[1];
  r  = circle.radius();
  const int maxIte = 2;
  for(int ite=0;ite<maxIte;++ite){
    Lpav circle2;
    circle2.clear();
    // MDC
    for(int i=0;i<list.length();++i){
      double R = sqrt((list[i]->wire()->xyPosition().x()-xc)*(list[i]->wire()->xyPosition().x()-xc)+
                      (list[i]->wire()->xyPosition().y()-yc)*(list[i]->wire()->xyPosition().y()-yc));
      if(R == 0.)continue;
      double U = 1./R;
      double dir = R > r ? -1. : 1.;
      double X = xc+(list[i]->wire()->xyPosition().x()-xc)*U*(R+dir*list[i]->hit()->drift());
      double Y = yc+(list[i]->wire()->xyPosition().y()-yc)*U*(R+dir*list[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();
    //std::cout << xc << ", " << yc << " : " << r << std::endl;
  }
  return true;
}

int
TBuilderCurl::stereoHit(double &xc, double &yc, double &r, double &q,
                        AList<TMLink> & list) const 
{
  if(list.length() == 0)return -1;

  HepPoint3D center(xc, yc, 0.);
  HepPoint3D tmp(-999., -999., 0.);
  for(unsigned i = 0, size = list.length(); i < size; ++i){
    TMDCWireHit &h = *const_cast<TMDCWireHit*>(list[i]->hit());
    HepVector3D X = 0.5*(h.wire()->forwardPosition() +
                      h.wire()->backwardPosition());
    HepVector3D x     = HepVector3D(X.x(), X.y(), 0.);
    HepVector3D w     = x - center;
    HepVector3D V     = h.wire()->direction();
    HepVector3D v     = HepVector3D(V.x(), V.y(), 0.);
    double   vmag2 = v.mag2();
    double   vmag  = sqrt(vmag2);
    //...temporary
    for(unsigned j = 0; j < 4; ++j)
      list[i]->arcZ(tmp,j);
    
    //...stereo?
    if (vmag == 0.) continue;
   
    double drift = h.drift(WireHitLeft);
    double R[2] = {r + drift, r - drift};
    double wv = w.dot(v);
    double d2[2];
    d2[0] = vmag2*R[0]*R[0] + (wv*wv - vmag2*w.mag2()); //...= v^2*(r^2 - w^2*sin()^2)...outer
    d2[1] = vmag2*R[1]*R[1] + (wv*wv - vmag2*w.mag2()); //...= v^2*(r^2 - w^2*sin()^2)...inner
    //...No crossing in R/Phi plane...
    if (d2[0] < 0. && d2[1] < 0.) continue;

    bool ok_inner(true);
    bool ok_outer(true);
    double d[2] = {-1., -1.};
    //...outer
    if(d2[0] >= 0.){
      d[0] = sqrt(d2[0]);
    }else{
      ok_outer = false;
    }
    if(d2[1] >= 0.){
      d[1] = sqrt(d2[1]);
    }else{
      ok_inner = false;
    }
    
    //...Cal. length and z to crossing points...
    double l[2][2];
    double z[2][2];
    //...outer
    if(ok_outer){
      l[0][0] = (- wv + d[0]) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
      l[1][0] = (- wv - d[0]) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
      z[0][0] = X.z() + l[0][0]*V.z();
      z[1][0] = X.z() + l[1][0]*V.z();
    }
    //...inner
    if(ok_inner){
      l[0][1] = (- wv + d[1]) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() + (r^2 - w^2*sin()^2)^0.5)/v
      l[1][1] = (- wv - d[1]) / vmag2; //...= (-wvcos()+d)/v/v = (-wcos() - (r^2 - w^2*sin()^2)^0.5)/v
      z[0][1] = X.z() + l[0][1]*V.z();
      z[1][1] = X.z() + l[1][1]*V.z();
    }
    
    //...Cal. xy position of crossing points...
    HepVector3D p[2][2];
    if(ok_outer){
      p[0][0] = x + l[0][0] * v;
      p[1][0] = x + l[1][0] * v;
      HepVector3D tmp_pc = p[0][0] - center;
      HepVector3D pc0 = HepVector3D(tmp_pc.x(), tmp_pc.y(), 0.);
      p[0][0] -= drift/pc0.mag()*pc0;
      tmp_pc = p[1][0] - center;
      HepVector3D pc1 = HepVector3D(tmp_pc.x(), tmp_pc.y(), 0.);
      p[1][0] -= drift/pc1.mag()*pc1;
    }
    if(ok_inner){
      p[0][1] = x + l[0][1] * v;
      p[1][1] = x + l[1][1] * v;
      HepVector3D tmp_pc = p[0][1] - center;
      HepVector3D pc0 = HepVector3D(tmp_pc.x(), tmp_pc.y(), 0.);
      p[0][1] += drift/pc0.mag()*pc0;
      tmp_pc = p[1][1] - center;
      HepVector3D pc1 = HepVector3D(tmp_pc.x(), tmp_pc.y(), 0.);
      p[1][1] += drift/pc1.mag()*pc1;
    }

    //...Check r-phi...
    bool ok_xy[2][2];
    if(ok_outer){
      ok_xy[0][0] = true;
      ok_xy[1][0] = true;
    }else{
      ok_xy[0][0] = false;
      ok_xy[1][0] = false;
    }
    if(ok_inner){
      ok_xy[0][1] = true;
      ok_xy[1][1] = true;
    }else{
      ok_xy[0][1] = false;
      ok_xy[1][1] = false;
    }
    if(ok_outer){
      if (q * (center.x() * p[0][0].y() - center.y() * p[0][0].x())  < 0.)
        ok_xy[0][0] = false;
      if (q * (center.x() * p[1][0].y() - center.y() * p[1][0].x())  < 0.)
        ok_xy[1][0] = false;
    }
    if(ok_inner){
      if (q * (center.x() * p[0][1].y() - center.y() * p[0][1].x())  < 0.)
        ok_xy[0][1] = false;
      if (q * (center.x() * p[1][1].y() - center.y() * p[1][1].x())  < 0.)
        ok_xy[1][1] = false;
    }
    if(!ok_inner && ok_outer && (!ok_xy[0][0]) && (!ok_xy[1][0])){
      continue;
    }
    if(ok_inner && !ok_outer && (!ok_xy[0][1]) && (!ok_xy[1][1])){
      continue;
    }

    //...Check z position...
//Liuqg060925, change temporary! these should be changed to bes3, reference is in TTrack::szPosition!!! 
    if(ok_xy[0][0]){
      if (z[0][0] < h.wire()->backwardPosition().z() || 
          z[0][0] > h.wire()->forwardPosition().z()) ok_xy[0][0] = false;
    }
    if(ok_xy[1][0]){
      if (z[1][0] < h.wire()->backwardPosition().z() || 
          z[1][0] > h.wire()->forwardPosition().z()) ok_xy[1][0] = false;
    }
    if(ok_xy[0][1]){
      if (z[0][1] < h.wire()->backwardPosition().z() || 
          z[0][1] > h.wire()->forwardPosition().z()) ok_xy[0][1] = false;
    }
    if(ok_xy[1][1]){
      if (z[1][1] < h.wire()->backwardPosition().z() || 
          z[1][1] > h.wire()->forwardPosition().z()) ok_xy[1][1] = false;
    }
    if ((!ok_xy[0][0]) && (!ok_xy[1][0]) &&
        (!ok_xy[0][1]) && (!ok_xy[1][1])){
      continue;
    }

    double cosdPhi, dPhi;
    //unsigned index = 0;
    unsigned indexL = 0;
    unsigned indexR = 0;
    //std::cout << "Stereo " << std::endl;
    // Q > 0 : Outer = L, Inner = R
    // Q < 0 : Outer = R, Inner = L
    if(ok_xy[0][0]){
      //...cal. arc length...
      cosdPhi = - center.dot((p[0][0] - center).unit()) / center.mag();
      if(fabs(cosdPhi) < 1.0){
        dPhi = acos(cosdPhi);
      }else if(cosdPhi >= 1.0){
        dPhi = 0.;
      }else{
        dPhi = M_PI;
      }
      //list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), index);
      //std::cout << r*dPhi << ", " << z[0][0] << std::endl;
      //++index;
      if(q > 0){
        //std::cout << "Outer L" << std::endl;
        if(indexL == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), 0);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), 2);
        ++indexL;
      }else{
        //std::cout << "Outer R" << std::endl;
        if(indexR == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), 1);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[0][0], 0.), 3);
        ++indexR;
      }
    }
    if(ok_xy[1][0]){
      //...cal. arc length...
      cosdPhi = - center.dot((p[1][0] - center).unit()) / center.mag();
      if(fabs(cosdPhi) < 1.0){
        dPhi = acos(cosdPhi);
      }else if(cosdPhi >= 1.0){
        dPhi = 0.;
      }else{
        dPhi = M_PI;
      }
      //list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), index);
      //std::cout << r*dPhi << ", " << z[1][0] << std::endl;
      //++index;
      if(q > 0){
        //std::cout << "Outer L" << std::endl;
        if(indexL == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), 0);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), 2);
        ++indexL;
      }else{
        //std::cout << "Outer R" << std::endl;
        if(indexR == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), 1);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[1][0], 0.), 3);
        ++indexR;
      }
    }
    if(ok_xy[0][1]){
      //...cal. arc length...
      cosdPhi = - center.dot((p[0][1] - center).unit()) / center.mag();
      if(fabs(cosdPhi) < 1.0){
        dPhi = acos(cosdPhi);
      }else if(cosdPhi >= 1.0){
        dPhi = 0.;
      }else{
        dPhi = M_PI;
      }
      //list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), index);
      //std::cout << r*dPhi << ", " << z[0][1] << std::endl;
      //++index;
      if(q > 0){
        //std::cout << "Inner R" << std::endl;
        if(indexR == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), 1);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), 3);
        ++indexR;
      }else{
        //std::cout << "Inner L" << std::endl;
        if(indexL == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), 0);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[0][1], 0.), 2);
        ++indexL;
      }
    }
    if(ok_xy[1][1]){
      //...cal. arc length...
      cosdPhi = - center.dot((p[1][1] - center).unit()) / center.mag();
      if(fabs(cosdPhi) < 1.0){
        dPhi = acos(cosdPhi);
      }else if(cosdPhi >= 1.0){
        dPhi = 0.;
      }else{
        dPhi = M_PI;
      }
      //list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), index);
      //std::cout << r*dPhi << ", " << z[1][1] << std::endl;
      //++index;
      if(q > 0){
        //std::cout << "Inner R" << std::endl;
        if(indexR == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), 1);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), 3);
        ++indexR;
      }else{
        //std::cout << "Inner L" << std::endl;
        if(indexL == 0)list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), 0);
        else list[i]->arcZ(HepPoint3D(r*dPhi, z[1][1], 0.), 2);
        ++indexL;
      }
    }
  }
  return 0;
}

void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
                      AList<TMLink> &alayer0,AList<TMLink> &alayer1,
                      unsigned ip) const {
  AList<TMLink> tmp = alayer0;
  tmp.append(alayer1);
  double xc,yc,r;
  if(fitWDD(xc,yc,r,tmp)){
    double q = track.charge();
    stereoHit(xc,yc,r,q,slayer);
  }
}

void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
                      AList<TMLink> &alayer0,AList<TMLink> &alayer1,
                      AList<TMLink> &alayer2,
                      unsigned ip) const {
  AList<TMLink> tmp = alayer0;
  tmp.append(alayer1);
  tmp.append(alayer2);
  double xc,yc,r;
  if(fitWDD(xc,yc,r,tmp)){
    double q = track.charge();
    stereoHit(xc,yc,r,q,slayer);
  }
}

void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
                      AList<TMLink> &alayer0,AList<TMLink> &alayer1,
                      AList<TMLink> &alayer2,AList<TMLink> &alayer3,
                      unsigned ip) const {
  AList<TMLink> tmp = alayer0;
  tmp.append(alayer1);
  tmp.append(alayer2);
  tmp.append(alayer3);
  double xc,yc,r;
  if(fitWDD(xc,yc,r,tmp)){
    double q = track.charge();
    stereoHit(xc,yc,r,q,slayer);
  }
}

void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
                      AList<TMLink> &alayer0,AList<TMLink> &alayer1,
                      AList<TMLink> &alayer2,AList<TMLink> &alayer3,
                      AList<TMLink> &alayer4,
                      unsigned ip) const {
  AList<TMLink> tmp = alayer0;
  tmp.append(alayer1);
  tmp.append(alayer2);
  tmp.append(alayer3);
  tmp.append(alayer4);
  double xc,yc,r;
  if(fitWDD(xc,yc,r,tmp)){
    double q = track.charge();
    stereoHit(xc,yc,r,q,slayer);
  }
}

/*void
TBuilderCurl::setArcZ(TTrack &track, AList<TMLink> &slayer,
                      AList<TMLink> &alayer0,AList<TMLink> &alayer1,
                      AList<TMLink> &alayer2,AList<TMLink> &alayer3,
                      AList<TMLink> &alayer4,AList<TMLink> &alayer5,
                      unsigned ip) const {
  AList<TMLink> tmp = alayer0;
  tmp.append(alayer1);
  tmp.append(alayer2);
  tmp.append(alayer3);
  tmp.append(alayer4);
  tmp.append(alayer5);
  double xc,yc,r;
  if(fitWDD(xc,yc,r,tmp)){
    double q = track.charge();
    stereoHit(xc,yc,r,q,slayer);
  }
}
*/
#undef DEBUG_CURL_DUMP
#undef DEBUG_CURL_GNUPLOT
#undef DEBUG_CURL_MC

// End === Stereo Finder For Curl Tracks : by jtanaka ===

---