MdcxFindTracks Class Reference

#include <MdcxFindTracks.h>

List of all members.

Public Member Functions
	MdcxFindTracks ()
	MdcxFindTracks (const HepAList< MdcxSeg > &f, int debug=0)
virtual	~MdcxFindTracks ()
const HepAList< MdcxFittedHel > &	GetMdcxTrklist ()
void	process (const HepAList< MdcxSeg > &f)
void	print (std::ostream &o)
void	plot () const
Protected Member Functions
void	KillList ()
void	beginmore ()
void	endmore ()
double	findz_sl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
double	findz_cyl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
double	dlen (int slay1, double p1, int slay2, double p2, double om)
MdcxFittedHel	drophits (MdcxFittedHel *f)
MdcxHel	TakeToOrigin (MdcxHel &)
void	resout (MdcxFittedHel *f)
void	printpar (std::ostream &o)
bool	testFromSameTrack (MdcxSeg *seg1, MdcxSeg *seg2)
Protected Attributes
HepAList< MdcxFittedHel >	MdcxTrklist
MdcxHel	mhel
float	probmin
float	curvmax
int	ngood
int	npure
int	nbad
int	nfail
int	nhitsmc
int	nhitsgd
int	m_debug
bool	m_doSag

---

Detailed Description

Definition at line 32 of file MdcxFindTracks.h.

---

Constructor & Destructor Documentation

MdcxFindTracks::MdcxFindTracks

(

)

Definition at line 110 of file MdcxFindTracks.cxx.

References curvmax, MdcxParameters::minTrkProb, and probmin.

                              { 
  probmin = MdcxParameters::minTrkProb;
  curvmax = 1000000.0;
}

MdcxFindTracks::MdcxFindTracks	(	const HepAList< MdcxSeg > &	f,
		int	debug = 0
	)

Definition at line 115 of file MdcxFindTracks.cxx.

References curvmax, m_debug, MdcxParameters::minTrkProb, probmin, and process().

                                                                         { 
  probmin = MdcxParameters::minTrkProb;
  curvmax = 1000000.0;
  m_debug = debug;
  process(segList); 
}

MdcxFindTracks::~MdcxFindTracks

(

)

[virtual]

Definition at line 122 of file MdcxFindTracks.cxx.

References KillList().

{ KillList(); }

---

Member Function Documentation

void MdcxFindTracks::beginmore

(

)

[protected]

Definition at line 839 of file MdcxFindTracks.cxx.

References nbad, nfail, ngood, nhitsgd, nhitsmc, npure, and printpar().

                              {
  ngood=0; nbad=0; nfail=0; npure=0; // zero scalars
  nhitsmc=0; nhitsgd=0;         // zero scalars
  printpar(cout);               // print parameters to screen
}//endof beginmore

double MdcxFindTracks::dlen	(	int	slay1,
		double	p1,
		int	slay2,
		double	p2,
		double	om
	)			[protected]

Definition at line 894 of file MdcxFindTracks.cxx.

References Constants::pi.

Referenced by process(), and TakeToOrigin().

                                                                                {
  //double dp = p2 - p1;
  double dp;//yzhang change 2011-10-17 
  if((slay1==2 || slay1==3 || slay1 == 4 || slay1 == 9 || slay1 ==10) ){ //Ax
    dp = p2 - p1;
  }else if((slay2 >= slay1)) {
    dp = p2 - p1;
  }else {
    dp = p1 - p2;
  }

  if (om < 0.0) {
    while (dp < -Constants::pi) dp += Constants::pi;
    while (dp > 0.0)   dp -= 1*Constants::pi;
  } else {
    while (dp < 0.0)   dp += 1*Constants::pi;
    while (dp > Constants::pi)  dp -= 1*Constants::pi;
  }
  //double dl = dp * r;
  double dl = dp/om;
  //if(4 == m_debug) {
  //  double r = 1./om;
  //  cout prt(p1) prt(p2) prt(dp) prt(om) prt(r) prt(dl) << endl;
  //}
  return dl;
}

MdcxFittedHel MdcxFindTracks::drophits

(

MdcxFittedHel *

f

)

[protected]

Definition at line 921 of file MdcxFindTracks.cxx.

References MdcxHel::Doca(), MdcxParameters::dropHitsSigma, g_eventNo, m_debug, m_segDropHitsDoca, m_segDropHitsDrift, m_segDropHitsEvtNo, m_segDropHitsLayer, m_segDropHitsMcTkId, m_segDropHitsPull, m_segDropHitsSigma, m_segDropHitsWire, m_xtupleDropHits, and MdcxFittedHel::XHitList().

Referenced by process().

                                                             {
  MdcxHel mdcxHel = (MdcxHel) *finehel;
  const HepAList<MdcxHit>& hitlist = finehel->XHitList();
  HepAList<MdcxHit> nwhitlist = hitlist;
  int nremove = 0;
  int kkk = 0;
  while(nwhitlist[kkk]) {
    double pull = nwhitlist[kkk]->pull(mdcxHel);
    int layer = nwhitlist[kkk]->Layer();

    if(m_xtupleDropHits){
      float t_doca = mdcxHel.Doca(*(nwhitlist[kkk]));
      float t_e =  nwhitlist[kkk]->e(t_doca);
      m_segDropHitsEvtNo = g_eventNo;
      m_segDropHitsLayer = layer;
      m_segDropHitsWire = nwhitlist[kkk]->WireNo();
      m_segDropHitsPull = pull;
      m_segDropHitsDoca = t_doca;
      m_segDropHitsSigma = t_e;
      m_segDropHitsDrift = nwhitlist[kkk]->d(mdcxHel); 
      m_segDropHitsMcTkId = nwhitlist[kkk]->getDigi()->getTrackIndex();
      m_xtupleDropHits->write();
    }

    if (fabs(pull) > MdcxParameters::dropHitsSigma[layer]){
      nwhitlist.remove(kkk);
      nremove++;
    } else {
      kkk++;
    }
  }
  if(m_debug==4) cout<< " MdcxFindTracks::drophits drop "<<nremove<<" hits "
    << " nhits of drop helix = "<<nwhitlist.length()<<endl;

  MdcxFittedHel newhel(nwhitlist,mdcxHel);
  return newhel;
}//endof drophits

void MdcxFindTracks::endmore

(

)

[protected]

Definition at line 845 of file MdcxFindTracks.cxx.

References print().

                            {
  print(cout);
}//endof endmore

double MdcxFindTracks::findz_cyl	(	int	i1,
		int	i2,
		const HepAList< MdcxSeg > &	slclus
	)			[protected]

Definition at line 869 of file MdcxFindTracks.cxx.

Referenced by process().

                                                                                    {
  double zint = -1000.0;
  double r = 1.0/fabs(segList[iAx]->Omega());
  double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0(); 
  double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
  double sx = segList[iStU]->Xline_slope();
  double sy = segList[iStU]->Yline_slope();
  double cx = xl_0 - segList[iAx]->Xc();
  double cy = yl_0 - segList[iAx]->Yc();
  double a = sx*sx + sy*sy;
  double b = 2.0*(sx*cx + sy*cy);
  //double c = (cx*cx + cy*cy - r*r);
  double bsq = b * b;
  double fac = 4.0 * a * (cx*cx + cy*cy - r*r);
  double ta = 2.0 * a;
  if (fac < bsq) {
    double sfac = sqrt(bsq - fac);
    double z1   = -(b - sfac)/ta;
    double z2   = -(b + sfac)/ta;
    zint = (fabs(z1) < fabs(z2)) ? z1 : z2;
  }
  // cout << " findz_cyl segs " << iAx << " " << iStU << " zint " << zint << endl;
  return zint;
}

double MdcxFindTracks::findz_sl	(	int	i1,
		int	i2,
		const HepAList< MdcxSeg > &	slclus
	)			[protected]

Definition at line 849 of file MdcxFindTracks.cxx.

References cos(), and sin().

Referenced by process().

                                                                                  {
  double Ap = -sin(segList[iAx]->Phi0_sl_approx());
  double Bp =  cos(segList[iAx]->Phi0_sl_approx());
  double xp = segList[iAx]->Xref() + Ap*segList[iAx]->D0_sl_approx();
  double yp = segList[iAx]->Yref() + Bp*segList[iAx]->D0_sl_approx();
  double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0(); 
  double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();

  const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
  double Cl = hitsStU[0]->wz();
  double ndotl = Bp*hitsStU[0]->wy() + Ap*hitsStU[0]->wx();
  /*
     std::cout<<"ndotl   "<<ndotl<<" "<<Bp<<" "<<hitsStU[0]->wy()<<" "<<
     Ap<<" "<<hitsStU[0]->wx()<<std::endl;
     */
  double zint = (ndotl != 0.0) ? ((Bp*(yp-yl)+Ap*(xp-xl))*Cl/ndotl) : -1000.;
  // cout << " findz_sl segs " << iAx << " " << iStU << " zint " << zint << endl;
  return zint;
}

const HepAList<MdcxFittedHel>& MdcxFindTracks::GetMdcxTrklist

(

)

[inline]

Definition at line 39 of file MdcxFindTracks.h.

References MdcxTrklist.

Referenced by MdcxTrackFinder::execute().

{return MdcxTrklist;}

void MdcxFindTracks::KillList

(

)

[inline, protected]

Definition at line 61 of file MdcxFindTracks.h.

References MdcxTrklist.

Referenced by ~MdcxFindTracks().

{HepAListDeleteAll(MdcxTrklist);}

void MdcxFindTracks::plot

(

)

const

Definition at line 830 of file MdcxFindTracks.cxx.

                               {
} // endof plot

void MdcxFindTracks::print

(

std::ostream &

o

)

Referenced by endmore().

void MdcxFindTracks::printpar

(

std::ostream &

o

)

[protected]

Referenced by beginmore().

void MdcxFindTracks::process

(

const HepAList< MdcxSeg > &

f

)

Definition at line 124 of file MdcxFindTracks.cxx.

References MdcxFittedHel::Chisq(), cos(), MdcxHel::D0(), dlen(), MdcxHel::Doca(), MdcxHel::Doca_Len(), drophits(), Constants::epsilon, MdcxFittedHel::Fail(), MdcxParameters::findTrkGroup, findz_cyl(), findz_sl(), MdcxFittedHel::FitPrint(), MdcxHel::flip(), g_dPhiAU, g_dPhiAV, g_eventNo, g_omegag, g_trklappend1, g_trklappend2, g_trklappend3, g_trklfirstProb, g_trklgood, g_trklhelix, g_trkllmk, g_trklproca, MdcxFittedHel::Grow(), MdcxParameters::helixFitSigma, genRecEmupikp::i, MdcxFittedHel::Layer(), MdcxParameters::layerSet2AddSeg, MdcxHel::Lmax(), m_addSegAddD0, m_addSegAddPhi, m_addSegAddPhi0, m_addSegAddPhiLay, m_addSegAddSl, m_addSegEvtNo, m_addSegLen, m_addSegPoca, m_addSegSame, m_addSegSeedD0, m_addSegSeedPhi, m_addSegSeedPhi0, m_addSegSeedPhiLay, m_addSegSeedSl, m_addSegSlayer, m_debug, m_segCombDLenAU, m_segCombDLenUV, m_segCombOmega, m_segCombPhiA, m_segCombPhiU, m_segCombPhiV, m_segCombSameAU, m_segCombSameUV, m_segCombSlA, m_segCombSlU, m_segCombSlV, m_xtupleAddSeg1, m_xtupleAddSeg2, m_xtupleSegComb, MdcxParameters::maxDlen, MdcxParameters::maxDp12, MdcxParameters::maxDp13, MdcxParameters::maxMdcZLen, MdcxParameters::maxProca, MdcxParameters::maxRcsInAddSeg, MdcxParameters::maxTrkOmega, Mdcxprobab(), MdcxTrklist, mhel, MdcxFittedHel::Nhits(), MdcxParameters::nSegCombo, Constants::nSuperLayer, MdcxHel::Omega(), MdcxHel::Phi0(), Constants::pi, MdcxHel::print(), MdcxFittedHel::Prob(), probmin, prt, MdcxFittedHel::Rcs(), MdcxFittedHel::ReFit(), sin(), SuperLayer(), TakeToOrigin(), testFromSameTrack(), Constants::twoPi, MdcxHel::X0(), MdcxHel::Xc(), MdcxFittedHel::XHitList(), MdcxHel::Y0(), and MdcxHel::Yc().

Referenced by MdcxFindTracks().

                                                             {
  static int   nfound  = 0;
  static float sumprob = 0.0; 

  MdcxHel thel(0., 0., 0., 0., 0.);
  mhel = thel;
  int nSeg = segList.length();
  for (int i = 0; i < nSeg ; i++) { segList[i]->SetUsedOnHel(0); }

  double xp, yp, xl1, yl1, rl1, xl2, yl2, rl2;
  double d0g, phi0g, omegag, z0g, tanlg;
  double d0g_sl, phi0g_sl, omegag_sl, z0g_sl, tanlg_sl;
  double zintAU, zintAV, zintAU_sl, zintAV_sl = 9999.;
  double rl1_sl, rl2_sl;
  double xc, yc;
  double sinPhi0g_sl, cosPhi0g_sl, xp_sl, yp_sl;
  double phiAx, phiStU, phiStV, xl1_sl, yl1_sl, xl2_sl, yl2_sl, xl1_0, yl1_0, xl2_0, yl2_0;
  double sx1, sy1, sx2, sy2;
  double x0g, y0g;  
  double fltAx, ztest, fltL_sl, ztest_sl;  
  int floop, trkCount= 0;

  //The main loop; first get the 3 SL combo we'll use to make a trial helix
  for (int iSegCombo = 0; iSegCombo < MdcxParameters::nSegCombo; iSegCombo++) {
    //FIXME optimize: group segList by superlayer ahead
    int axlay = MdcxParameters::findTrkGroup[iSegCombo][0];
    int stUlay = MdcxParameters::findTrkGroup[iSegCombo][1];
    int stVlay = MdcxParameters::findTrkGroup[iSegCombo][2];

    if(4 == m_debug) std::cout <<" Test combo: ("<< axlay<<","<<stUlay<<","<<stVlay <<")"<< std::endl; 

    float max_dPhiAU = MdcxParameters::maxDp12[iSegCombo];
    float max_dPhiAV = MdcxParameters::maxDp13[iSegCombo];

    //---------------------------------------------------------------------
    //  Axial SL segment furnishes initial guess for d0, phi0, omega
    //---------------------------------------------------------------------
    for (int iAx = 0; iAx < nSeg; iAx++) {
      bool b_useGood_stU = true; //yzhang add 2009-10-21 17:20:45
      bool b_useGood_stV = true; //yzhang add 2009-10-21 17:20:45
      if ((segList[iAx]->GetUsedOnHel() != 0) || (segList[iAx]->Fail() != 0) 
          || (segList[iAx]->SuperLayer(0) != axlay) || (segList[iAx]->Origin() != -1) ) continue;

      if (4 == m_debug){
        std::cout<< "---------1.Ax seg------ No."<<iAx<< std::endl;
        segList[iAx]->printSegAll();
      }

      //Skip low pt track
      omegag = segList[iAx]->Omega();
      if (fabs(omegag) > MdcxParameters::maxTrkOmega) {//FIXME 0.2
        if (4 == m_debug) std::cout<< "SKIP by maxTrkOmega "<<fabs(omegag) << std::endl;
        continue;  
      }
      if(4 == m_debug) std::cout << " iAx omegag = " << omegag << std::endl;
      if(g_omegag) g_omegag->fill(omegag);

      const HepAList<MdcxHit>& hitsAx = segList[iAx]->XHitList();

      //Interpret segment as straight line also
      d0g_sl = segList[iAx]->D0_sl_approx();
      phi0g_sl = segList[iAx]->Phi0_sl_approx(); 
      omegag_sl = 0.0;

      //save c-tors; do it locally
      sinPhi0g_sl = -sin(phi0g_sl);
      cosPhi0g_sl = cos(phi0g_sl);
      xp_sl = segList[iAx]->Xref() + sinPhi0g_sl*d0g_sl;
      yp_sl = segList[iAx]->Yref() + cosPhi0g_sl*d0g_sl;
      d0g_sl = yp_sl*cosPhi0g_sl + xp_sl*sinPhi0g_sl;
      fltL_sl = segList[iAx]->Xref()*cosPhi0g_sl - segList[iAx]->Yref()*sinPhi0g_sl; 
      if(4 == m_debug) {
        std::cout << "--Ax :" prt(d0g_sl) prt(phi0g_sl) prt(omegag_sl) prt(sinPhi0g_sl) prt(xp_sl) prt(yp_sl) prt(fltL_sl) << std::endl; 
      }

      if (fabs(omegag) > Constants::epsilon) {
        MdcxHel ohel = TakeToOrigin(*segList[iAx]);
        omegag = ohel.Omega();
        phi0g = ohel.Phi0();
        d0g = ohel.D0(); 
        xc = ohel.Xc();
        yc = ohel.Yc();
        x0g = ohel.X0();
        y0g = ohel.Y0();
        phiAx = atan2(y0g-yc, x0g-xc);
        fltAx = dlen(-2, phi0g,-1, segList[iAx]->Phi0(), omegag); 
        //fltAx = dlen(phi0g,segList[iAx]->Phi0(), omegag); //yzhang TEMP 2011-10-17 
        if(4 == m_debug) {
          std::cout <<"--Ax :"<< " ohel ";
          ohel.print();
        }
      }

      //---------------------------------------------------------------------
      //First stereo SL segement and axial segment furnish one z intersection
      //---------------------------------------------------------------------
      //for (int iStU = 0; iStU < nSeg; iStU++ ) //yzhang delete 2009-10-21 18:13:50
      for (int iStU = -1; true; ) {//yzhang add
        if(!b_useGood_stU && (iStU >= (nSeg-1))) break;
        if(b_useGood_stU && (iStU >= (nSeg-1))){
          iStU = -1;
          b_useGood_stU = false;
        }
        iStU++;
        if ((segList[iAx]->GetUsedOnHel() != 0) 
            || (segList[iStU]->GetUsedOnHel() != 0) 
            || (segList[iStU]->Fail() != 0) 
            || (segList[iStU]->SuperLayer(0) != stUlay) 
            || (segList[iStU]->Origin() != -1) 
            //yzhang add  2009-10-21 17:21:55
            || (b_useGood_stU && (segList[iStU]->Nhits()< 4))
           ){ continue;}

        if (4 == m_debug){
          std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl; 
          std::cout<< "---------2.St U seg ------No."<<iStU<< std::endl;//yzhang debug
          segList[iStU]->printSeg();
          std::cout<<" omega of seg iStU = "<<segList[iStU]->Omega()<<std::endl;
        }

        //Test dPhi between Ax slayer and StU slayer 
        const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
        double dPhiAU = fabs(hitsAx[0]->phiMid() - hitsStU[0]->phiMid());
        if (dPhiAU > Constants::pi) dPhiAU = Constants::twoPi - dPhiAU;
        if(4 == m_debug) {
          if (dPhiAU > max_dPhiAU){
            cout << "SKIP by dPhiAU " <<dPhiAU <<" > "<< max_dPhiAU << endl;
          }else{
            cout << "KEEP by dPhiAU " <<dPhiAU<< " <= " << max_dPhiAU << endl;
          }
          std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
            <<" StU "<<hitsStU[0]->phiMid() << std::endl; 
          std::cout<<  std::endl;
        }
        if (g_dPhiAU) g_dPhiAU->fill(dPhiAU,segList[iStU]->SuperLayer());

        if (dPhiAU > max_dPhiAU){ continue; }

        //Calculate z by Ax and StU segs
        xl1_0 = segList[iStU]->Xline_bbrrf();
        yl1_0 = segList[iStU]->Yline_bbrrf();
        sx1   = segList[iStU]->Xline_slope();
        sy1   = segList[iStU]->Yline_slope();
        if (4 == m_debug) std::cout prt(xl1_0) prt(yl1_0) prt(sx1) prt(sy1) prt(omegag) << std::endl;
        if (fabs(omegag) >Constants::epsilon) {
          zintAU = findz_cyl(iAx, iStU, segList);
          xl1 = xl1_0 + sx1*zintAU;
          yl1 = yl1_0 + sy1*zintAU;
          rl1 = sqrt(xl1*xl1 + yl1*yl1);
          phiStU = atan2(yl1-yc, xl1-xc);
          if (4 == m_debug) cout prt(zintAU) prt(xl1) prt(yl1) prt(rl1) prt(phiStU) <<endl;
        } else {
          zintAU = -9999.;
        }

        //Calculate z_sl by Ax and StU segs
        zintAU_sl = findz_sl(iAx, iStU, segList); 
        xl1_sl = xl1_0 + sx1*zintAU_sl;
        yl1_sl = yl1_0 + sy1*zintAU_sl;
        rl1_sl = sqrt(xl1_sl*xl1_sl + yl1_sl*yl1_sl);
        if (4 == m_debug) cout prt(zintAU_sl) prt(xl1_sl) prt(yl1_sl) prt(rl1_sl) <<endl;

        //Cut by z and z_sl
        if ( (zintAU < -MdcxParameters::maxMdcZLen) &&
            (zintAU_sl < -MdcxParameters::maxMdcZLen) ) {
          if (4 == m_debug) std::cout << " SKIP by zintAU < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
          continue;
        }
        if ( (zintAU > MdcxParameters::maxMdcZLen) &&
            (zintAU_sl > MdcxParameters::maxMdcZLen) ) {
          if (4 == m_debug) std::cout << " SKIP by zintAU > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
          continue;
        }

        //---------------------------------------------------------------------
        //Next stereo SL segement and axial segment furnish other z intersection
        //---------------------------------------------------------------------
        for (int iStV = -1; true; ) {
          if(!b_useGood_stV && (iStV >= (nSeg-1))) {
            if (4 == m_debug) std::cout << iStV <<" no segments in V slayer "<< std::endl; 
            break;
          }
          if(b_useGood_stV && (iStV >= (nSeg-1))){
            iStV = -1;
            b_useGood_stV = false;
          }
          iStV++;
          if ((segList[iStV]->GetUsedOnHel()!=0) 
              || (segList[iStU]->GetUsedOnHel()!=0) 
              || (segList[iAx]->GetUsedOnHel()!=0) 
              || (segList[iStV]->Fail() != 0) 
              || (segList[iStV]->SuperLayer(0) != stVlay) 
              || (segList[iStV]->Origin() != -1)  
              //yzhang add  2009-10-21 17:21:55
              || (b_useGood_stV && (segList[iStV]->Nhits()< 4))
             ){ continue; }

          if (4 == m_debug){
            std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl; 
            std::cout<< "---------3.St V seg ------No."<<iStV<< std::endl;//yzhang debug
            segList[iStV]->printSeg();
            std::cout<<" omega of seg iStV = "<<segList[iStV]->Omega()<<std::endl;
          }

          //Calculate dPhi between Ax and StV segs
          const HepAList<MdcxHit>& hitsStV = segList[iStV]->XHitList();
          double dPhiAV = fabs(hitsAx[0]->phiMid() - hitsStV[0]->phiMid());
          if (dPhiAV > Constants::pi) dPhiAV = Constants::twoPi - dPhiAV;

          if(4 == m_debug) {
            if (dPhiAV > max_dPhiAV){
              cout << "SKIP by dPhiAV " <<dPhiAV <<" > "<< max_dPhiAV << endl;
            }else{
              cout << "KEEP by dPhiAV " <<dPhiAV<< " <= " << max_dPhiAV << endl;
            }
            std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
              <<" StV "<<hitsStV[0]->phiMid() << std::endl; 
            std::cout<<  std::endl;
          }
          if (g_dPhiAV) g_dPhiAV->fill(dPhiAV,segList[iStV]->SuperLayer());
          if (dPhiAV > max_dPhiAV) { continue; }


          //Calculate z by Ax and StV segs
          xl2_0 = segList[iStV]->Xline_bbrrf();
          yl2_0 = segList[iStV]->Yline_bbrrf();
          sx2   = segList[iStV]->Xline_slope();
          sy2   = segList[iStV]->Yline_slope();
          if (fabs(omegag) > Constants::epsilon) {
            zintAV = findz_cyl(iAx, iStV, segList);
            xl2 = xl2_0 + sx2*zintAV;
            yl2 = yl2_0 + sy2*zintAV;
            rl2 = sqrt(xl2*xl2 + yl2*yl2);
            if (4 == m_debug){ 
              segList[iAx]->printSeg();
              segList[iStV]->printSeg();
              cout << "zintAV " << zintAV << endl;
            }
            phiStV = atan2(yl2-yc, xl2-xc);
          } else {
            zintAV = -9999.;
          }

          //Calculate z by Ax and StV segs
          zintAV_sl = findz_sl(iAx, iStV, segList);
          xl2_sl = xl2_0 + sx2*zintAV_sl;
          yl2_sl = yl2_0 + sy2*zintAV_sl; 
          rl2_sl = sqrt(xl2_sl*xl2_sl + yl2_sl*yl2_sl);

          //Cut by z of AV segs
          if ( (zintAV < -MdcxParameters::maxMdcZLen) &&
              (zintAV_sl < -MdcxParameters::maxMdcZLen) ) {
            if (4 == m_debug) std::cout << " SKIP by zintAV < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
            continue;
          }
          if ( (zintAV > MdcxParameters::maxMdcZLen) &&
              (zintAV_sl > MdcxParameters::maxMdcZLen) ) {
            if (4 == m_debug) std::cout << " SKIP by zintAV > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
            continue;
          }

          //We now have enough info to form and fit a guess helix (inside iAx-iStU-iStV loops)
          if (4 == m_debug) {
            //cout << "KEEP layer[klm]: [" << axlay <<","<< stUlay ","<< stVlay << "] " prt(d0g) prt(phi0g) prt(omegag) zrt(zintAU) prt(zintAV) << endl; 
            segList[iAx]->printSeg();
            segList[iStU]->printSeg();
            segList[iStV]->printSeg();
          }
          MdcxFittedHel fithel;
          double helixFitSigma = MdcxParameters::helixFitSigma;
          double first_prob = 0.0; 
          HepAList<MdcxHit> hitlist;
          hitlist.append(hitsAx);
          hitlist.append(hitsStU);
          hitlist.append(hitsStV);
          if(g_trkllmk) for (int i=0; i<hitlist.length(); i++ ){ g_trkllmk->fill(hitlist[i]->Layer()); }


          //---------------------------------------------------------------------
          //Try fitting the trial helix from axial SL as circle first
          //---------------------------------------------------------------------
          floop = 1;
          while (floop) {
            if(4 == m_debug) std::cout<< "---------4.Ax circle fit------"<< std::endl;//yzhang debug
            if (fabs(omegag) < Constants::epsilon) {
              if(4 == m_debug) std::cout<<"SKIP by omegag==0   "<<std::endl;
              break;
            }
            if ( (zintAU < -MdcxParameters::maxMdcZLen) ||
                (zintAU > MdcxParameters::maxMdcZLen) ){
              if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
              break; 
            }
            if ( (zintAV < -MdcxParameters::maxMdcZLen) ||
                (zintAV > MdcxParameters::maxMdcZLen) ){
              if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
              break; 
            }


            //Calculate flight length of AUV segs
            if(4 == m_debug) cout<< "dlen calc. slay U "<<segList[iStU]->SuperLayer()<< " slay V "<<segList[iStV]->SuperLayer()<<endl;
            double fltLenUV = dlen(segList[iStU]->SuperLayer(), phiStU,segList[iStV]->SuperLayer(), phiStV, omegag);//yzhang change TEMP 2011-10-17 
            //double fltLenUV = dlen(phiStU, phiStV, omegag);
            if(4 == m_debug){
              cout prt(fltLenUV) prt(phiStU) prt(phiStV) prt(omegag)<< std::endl;
            }
            if (fltLenUV > MdcxParameters::maxDlen) {//15. to 150.
              if(4 == m_debug) std::cout<<"SKIP by fltLenUV > "<<MdcxParameters::maxDlen<<std::endl;
              break;   //FIXME
            }
            tanlg = (zintAV - zintAU)/fltLenUV; 
            if(4 == m_debug) cout<< "dlen calc. slay A "<<segList[iAx]->SuperLayer()<< " slay U "<<segList[iStU]->SuperLayer()<<endl;
            double fltLenAU = dlen(segList[iAx]->SuperLayer(), phiAx, segList[iStU]->SuperLayer(), phiStU, omegag);//yzhang change TEMP 2011-10-17 

            if(m_xtupleSegComb){
              m_segCombOmega = omegag;
              m_segCombSameAU = testFromSameTrack(segList[iAx],segList[iStU]);
              m_segCombSameUV = testFromSameTrack(segList[iStU],segList[iStV]);
              m_segCombSlA= segList[iAx]->SuperLayer();
              m_segCombSlU= segList[iStU]->SuperLayer();
              m_segCombSlV= segList[iStV]->SuperLayer();
              m_segCombPhiA = phiAx;
              m_segCombPhiU = phiStU;
              m_segCombPhiV = phiStV;
              m_segCombDLenUV = fltLenUV;
              m_segCombDLenAU = fltLenAU;
              m_xtupleSegComb->write();
            }
            //double fltLenAU = dlen(phiAx, phiStU, omegag);
            z0g = zintAU - fltLenAU*tanlg;
            if (4 == m_debug) cout prt(z0g) prt(tanlg) prt(fltLenAU) prt(zintAU) prt(zintAV)<< endl;

            //Cut by z0g
            if ((z0g < -MdcxParameters::maxMdcZLen) || (z0g > MdcxParameters::maxMdcZLen)) {
              if (4 == m_debug) std::cout<<"SKIP by z0g out of range "<<z0g<<">"<<MdcxParameters::maxMdcZLen<<std::endl;
              break;
            }
            ztest = z0g + fltAx*tanlg; 
            if (4 == m_debug) std::cout prt(ztest) prt(fltAx)<<std::endl;

            //---Helix fitting
            //MdcxHel ghel(d0g,phi0g,omegag,z0g,tanlg); // ghel.print();
            MdcxHel ghel(segList[iAx]->D0(), segList[iAx]->Phi0(), segList[iAx]->Omega(),
                ztest, tanlg, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref()); 
            //zoujh?: ghel.SetTurnFlag(1); // ghel.print();
            MdcxFittedHel firstfit(hitlist,ghel,helixFitSigma);  
            first_prob = firstfit.Prob(); 
            if (4 == m_debug) {
              std::cout<<" after first fit:  "<<std::endl;
              firstfit.FitPrint();
            }
            if (first_prob >= probmin) {
              MdcxHel helixOrigin= TakeToOrigin(firstfit);
              MdcxFittedHel nextfit(hitlist,helixOrigin,helixFitSigma);
              first_prob = nextfit.Prob();
              fithel=nextfit;
              if (g_trklgood) for (int i=0; i<nextfit.Nhits(); i++){ g_trklgood->fill(nextfit.Layer(i));}
            }
            if ( g_trklfirstProb) g_trklfirstProb->fill(first_prob);
            if (4 == m_debug) cout << " prob after helix fitting = " << first_prob << endl;     
            floop = 0;
          }//end fitting trial helix

          //---------------------------------------------------------------------
          // Try fitting straight line trial helix if fit from axial-SL-as-circle trial
          //  if helix not good enough
          //---------------------------------------------------------------------
          floop = 1;
          while (floop) {
            if (4 == m_debug) std::cout<< "---------4.2 straight line fit------"<< std::endl;//yzhang debug
            if (4 == m_debug) cout << " zintAU_sl " << zintAU_sl << " zintAV_sl " << zintAV_sl << endl; 
            if (first_prob > probmin) {
              if (4 == m_debug) cout << " helix fit good , exit straight line fit." << endl;     
              break;
            }

            if ( (zintAU_sl < -MdcxParameters::maxMdcZLen) ||
                (zintAU_sl > MdcxParameters::maxMdcZLen) ) break;
            if ( (zintAV_sl < -MdcxParameters::maxMdcZLen) ||
                (zintAV_sl > MdcxParameters::maxMdcZLen) ) break;

            double dx = xl2_sl - xl1_sl;
            double dy = yl2_sl - yl1_sl;
            double dl = sqrt(dx*dx + dy*dy);
            tanlg_sl = (zintAV_sl - zintAU_sl)/dl;
            dx = xl1_sl + d0g_sl*sin(phi0g_sl);
            dy = yl1_sl - d0g_sl*cos(phi0g_sl);
            dl = sqrt(dx*dx + dy*dy);
            z0g_sl = zintAU_sl - dl*tanlg_sl;
            if (4 == m_debug) cout prt(d0g_sl) prt(phi0g_sl) prt(z0g_sl) prt(tanlg_sl)<< endl;

            if ((z0g_sl < -MdcxParameters::maxMdcZLen) || (z0g_sl > MdcxParameters::maxMdcZLen)){
              if(4 == m_debug) std::cout << "SKIP by z0g_sl:"<<z0g_sl<<" > "<<MdcxParameters::maxMdcZLen << std::endl; 
              break;
            }
            //MdcxHel ghel(d0g_sl,phi0g_sl,omegag_sl,z0g_sl,tanlg_sl); // ghel.print();
            ztest_sl = z0g_sl + fltL_sl*tanlg_sl;
            if (4 == m_debug) cout prt(ztest_sl) << endl;

            //Helix fitting
            MdcxHel ghel(segList[iAx]->D0_sl_approx(), segList[iAx]->Phi0_sl_approx(), omegag_sl,
                ztest_sl, tanlg_sl, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref());
            //zoujh?: ghel.SetTurnFlag(1); // ghel.print();
            MdcxFittedHel firstfit(hitlist, ghel, helixFitSigma);
            first_prob = firstfit.Prob();
            //if (first_prob >= probmin)fithel=firstfit;
            if(4 == m_debug){
              ghel.print();
              std::cout <<"first_prob "<<first_prob << std::endl;  
            }

            if (first_prob >= probmin) {
              MdcxHel helixOrigin = TakeToOrigin(firstfit);
              MdcxFittedHel nextfit(hitlist, helixOrigin, helixFitSigma);
              first_prob = nextfit.Prob();
              fithel = nextfit;
              if (4 == m_debug) {
                cout << " three seg sl nextfit" <<  endl;
                nextfit.FitPrint();
              }
            }
            floop = 0;
          }

          //-----------------------------------------------------------------
          //Got a good trial helix fit; now try to add other segments onto it
          //-----------------------------------------------------------------
          floop = 1;
          while (floop) {
            floop = 0;
            if(4 == m_debug)std::cout<< "---------5. try add seg to helix------"<< std::endl;
            if (first_prob < probmin) {
              if(4 == m_debug) std::cout<< "prob"<<first_prob<< " "
                <<probmin<<" fit NOT good, exit add segs "<< std::endl;
              break;
            }
            float bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
            int   bndof = fithel.Nhits() - 5;
            float bprob = Mdcxprobab(bndof, bchisq);
            trkCount++; 
            segList[iAx]->SetUsedOnHel(trkCount);
            segList[iStU]->SetUsedOnHel(trkCount);
            segList[iStV]->SetUsedOnHel(trkCount); 

            if (4 == m_debug){
              cout << "in add seg to trail helix, klm seg:" << endl;  
              segList[iAx]->printSeg();                  
              segList[iStU]->printSeg();                  
              segList[iStV]->printSeg();                  
            }

            //Loop superlayers
            for (int iSlay = 0; iSlay < Constants::nSuperLayer; iSlay++) {
              int ilay = MdcxParameters::layerSet2AddSeg[iSegCombo][iSlay]; 
              for (int iSeg = 0; iSeg < nSeg; iSeg++) {
                if ((segList[iSeg]->SuperLayer(0) != ilay) 
                    ||(segList[iSeg]->GetUsedOnHel() != 0)
                    || (segList[iSeg]->Fail() != 0) 
                    || (segList[iSeg]->Origin() != -1)) continue;
                if(4 == m_debug) {
                  std::cout<< endl<< "== add seg "<< iSeg<<" ";
                  segList[iSeg]->printSeg();
                }

                //Cut by phi_diff = (Phi_Ax - Phi_Add)
                float phiAxSeg = segList[iAx]->XHitList()[0]->phiMid();
                float phiAddSeg = segList[iSeg]->XHitList()[0]->phiMid();
                double phiKNSegDiff = fabs(phiAxSeg - phiAddSeg);
                //yzhang add 2011-10-11 
                //double phiKNSegDiff = (phiAxSeg - phiAddSeg);
                //if(phiKNSegDiff<-Constants::pi) phiKNSegDiff+=2*Constants::pi;
                //if(phiKNSegDiff>Constants::pi) phiKNSegDiff-=2*Constants::pi;
                int t_same=-999;
                if(m_xtupleAddSeg1){
                  //std::cout<< __FILE__ << " testFromSameTrack  " << fromSameTrack<< "   "<<std::endl;
                  t_same = testFromSameTrack(segList[iAx],segList[iSeg]);
                  m_addSegSame = t_same;
                  m_addSegSeedSl = segList[iAx]->SuperLayer();
                  m_addSegSeedPhi = segList[iAx]->XHitList()[0]->phiMid();
                  m_addSegSeedPhiLay = segList[iAx]->XHitList()[0]->Layer();
                  m_addSegSeedPhi0 = segList[iAx]->Phi0_sl_approx();
                  m_addSegSeedD0 = segList[iAx]->D0_sl_approx();
                  m_addSegAddSl = segList[iSeg]->SuperLayer();
                  m_addSegAddPhi = segList[iSeg]->XHitList()[0]->phiMid();
                  m_addSegAddPhiLay = segList[iSeg]->XHitList()[0]->Layer();
                  m_addSegAddPhi0 = segList[iSeg]->Phi0_sl_approx();
                  m_addSegAddD0 = segList[iSeg]->D0_sl_approx();
                  m_xtupleAddSeg1->write();
                }
                //yzhang change 2011-10-11 
                if (phiKNSegDiff > 0.5*Constants::pi && phiKNSegDiff < 1.5*Constants::pi ) 
                  //if (fabs(phiKNSegDiff) > 0.7*Constants::pi ) 
                {//yzhang TEMP
                  if(4 == m_debug) std::cout<< " SKIP by phi diff,same "<< t_same <<" Ax "<<
                    phiAxSeg<<" iSeg "<<phiAddSeg<<" diff "<<phiKNSegDiff << std::endl;//yzhang debug 
                  continue;  // zoujh FIXME
                }else{
                  if(4 == m_debug) std::cout<< " phi diff OK, Ax "<<
                    phiAxSeg<<" added "<<phiAddSeg<<" diff="<<phiKNSegDiff << std::endl;//yzhang debug 
                }

                //Calculate doca
                xp = segList[iSeg]->Xref() - segList[iSeg]->SinPhi0()*segList[iSeg]->D0();
                yp = segList[iSeg]->Yref() + segList[iSeg]->CosPhi0()*segList[iSeg]->D0();
                const HepAList<MdcxHit>& hitsSegAdd= segList[iSeg]->XHitList();
                double proca = fithel.Doca( hitsSegAdd[0]->wx(), hitsSegAdd[0]->wy(),
                    hitsSegAdd[0]->wz(), xp, yp );
                if (g_trklproca) g_trklproca->fill(proca);

                //Segment passes loose doca cut; does it belong on this track?
                if(m_xtupleAddSeg2){
                  m_addSegEvtNo = g_eventNo;
                  m_addSegPoca = proca;
                  m_addSegSlayer = iSlay;
                  m_addSegLen = fithel.Doca_Len();
                  m_xtupleAddSeg2->write();
                }

                if (!((fabs(proca) < MdcxParameters::maxProca)&&(fithel.Doca_Len()<fithel.Lmax())) ){
                  //if (g_trklprocaSl) g_trklprocaSl->fill(segList[iSeg]->SuperLayer(0));
                  if(4 == m_debug){
                    std::cout<< " SKIP by fabs(proca):"<< fabs(proca)<< "<"<<MdcxParameters::maxProca <<" or Doca_Len:"<<fithel.Doca_Len() <<"<"<<fithel.Lmax()<< std::endl;//yzhang debug 
                  }
                }else{
                  if(4 == m_debug){
                    std::cout<< " proca & len OK, |proca| "<< fabs(proca)<< "<"<<MdcxParameters::maxProca <<" && Doca_Len:"<<fithel.Doca_Len() <<"<"<<fithel.Lmax()<< std::endl;//yzhang debug 
                  }
                }
                if(fithel.Doca_Len()<0) continue;//yzhang add TEMP 2011-10-11 

                if ((fabs(proca)<MdcxParameters::maxProca)&&(fithel.Doca_Len()<fithel.Lmax())) {
                  MdcxFittedHel newhel;
                  newhel.Grow(fithel, hitsSegAdd);
                  if (newhel.Nhits() != hitlist.length()) {
                    if (4 == m_debug){
                      cout<<" rcs "<<newhel.Rcs()<<"<=?"<<MdcxParameters::maxRcsInAddSeg<< endl;
                    }
                    //yzhang change 2009-10-21 FIXME
                    //if (newhel.Rcs() > MdcxParameters::maxRcsInAddSeg)
                    if (newhel.Prob() < probmin){
                      if(4 == m_debug){ 
                        cout << " prob " << newhel.Prob() << "<"<<probmin<<", ReFit" << endl;  
                      }
                      newhel.ReFit(); //FIXME
                    }
                    if(4 == m_debug) {
                      cout<<" refit prob "<<newhel.Prob()<<"<"<<probmin<<" Fail? "<<newhel.Fail()<< endl;
                    }
                    //yzhang change 2009-10-21 FIXME
                    if ( (newhel.Prob() >= probmin) && (newhel.Fail() == 0) ) {
                      //if ( (newhel.Rcs() <= MdcxParameters::maxRcsInAddSeg) && (newhel.Fail() == 0) ) 
                      fithel = newhel;
                      hitlist = newhel.XHitList();
                      segList[iSeg]->SetUsedOnHel(trkCount);
                      bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
                      bndof = fithel.Nhits() - 5;
                      float tprob = Mdcxprobab(bndof, bchisq); 
                      if (tprob > bprob) bprob = tprob;
                      if (4 == m_debug) {
                        cout << " segment ADDED, prob=" << newhel.Prob() << endl;
                      }
                    }   //*end trial track still good so keep new seg 
                    else{
                      if(4 == m_debug)std::cout<<" fit FAILED"<<std::endl;
                    }
                  } else { //*end new hits so refit
                    segList[iSeg]->SetUsedOnHel(trkCount);
                    if (4 == m_debug) cout << " segment ADDED, but no new hits" << endl;
                  }//*end dont need to refit if no new hits
                }//*end fabs(proca) < 0.010
              }//*end iSeg (loop over segs in a superlayer)
            }//*end iSlay (loop over superlayers)
            if (((fithel.Prob() < 0.90) && (fithel.Nhits() == 12)) || (fithel.Fail() != 0)) {
              for (int i = 0; i < nSeg; i++) { if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0); }
              trkCount--;
              break; //*jump to end floop first_prob >= probmin
            }

            if(4 == m_debug){
              std::cout<< "Track after add segs" << std::endl;//yzhang debug
              fithel.FitPrint();
            }

            //See if helix better moving backwards
            if(4 == m_debug) std::cout<< "---------6. flip------"<< std::endl;//yzhang debug
            fithel.flip();

            //See if this track shares lots of hits with another saved track
            int kki = 0;
            int notduplicate = 1; 
            while(MdcxTrklist[kki]) {
              if (4 == m_debug)std::cout<< "---------7. test saved track------"<< std::endl;//yzhang debug
              const HepAList<MdcxHit>& junk=MdcxTrklist[kki]->XHitList();
              int overl = 0; 
              int kkj = 0;
              while (junk[kkj]) {
                int kkl = 0;
                while (hitlist[kkl]) {
                  if (hitlist[kkl] == junk[kkj]) overl++; 
                  kkl++;
                }
                kkj++;
              }
              if (4 == m_debug) cout << "overlap with track # " << kki << " is "
                << junk.length() << " hitList " << hitlist.length() << " overlap " << overl << endl;
              if ( (hitlist.length()/4.) <= overl) notduplicate = 0;
              //yzhang 2009-10-21 16:30:53
              //if ( hitlist.length() <= (overl+2) ) notduplicate = 0;
              kki++;
            }

            if(g_trklhelix) for (int iHit=0; iHit<fithel.Nhits(); iHit++){ g_trklhelix->fill(fithel.Layer(iHit)); }

            //If not a duplicate, is it worth saving?
            if (notduplicate) {
              if(4 == m_debug) std::cout<< "---------8. test worth saving?------"<< std::endl;//yzhang debug
              MdcxFittedHel *finehel = new MdcxFittedHel(hitlist, fithel);

              if(4 == m_debug) {
                std::cout<<__FILE__
                  <<" finehel Prob>0.001 "<<finehel->Prob()
                  <<" nhits "<<finehel->Nhits() <<">=15? "
                  <<" bprob "<< bprob 
                  <<" >=?probmin "<< probmin
                  <<" Fail==0? "<< finehel->Fail()
                  <<  std::endl;
                finehel->FitPrint();
              }
              //cut nhits from 15->9 yzhang 2009-10-21 
              if (( (finehel->Prob()>0.001) || (finehel->Nhits()>=15) || (bprob>probmin) 
                    //|| (bchisq/bndof > MdcxParameters::maxRcsInAddSeg) 
                  ) && (finehel->Fail() == 0) ) {
                nfound++;
                sumprob += finehel->Prob();
                //if ( hitlist.length() == 16 )resout(finehel);

                //Begin hit dropping section
                //Set nodrop != 0 if want to prohibit dropping bad hits 
                int nodrop = 0; 
                if ((finehel->Prob() > 0.05) || nodrop ) {
                  MdcxTrklist.append(finehel);
                  //yzhang hist
                  for (int iHit=0; iHit<finehel->Nhits(); iHit++){
                    if(g_trklappend1) g_trklappend1->fill(finehel->Layer(iHit));//yzhang hist
                  }
                  //zhangy hist
                } else {
                  MdcxFittedHel* drophel = new MdcxFittedHel( drophits(finehel) );
                  float ratdrop = float(drophel->Nhits()) / float(finehel->Nhits());
                  if (4 == m_debug) cout << "APPEND track " << trkCount << ", drops " 
                    << finehel->Nhits()-drophel->Nhits() 
                      << " helixnhit "<<finehel->Nhits()<< " drophit "<<drophel->Nhits()
                      << " hits, drop rate="<<ratdrop <<">?"<<0.5
                      << " prob "<<drophel->Prob() <<" >?"<<finehel->Prob()
                      << " fail==0? "<<drophel->Fail()
                      << endl;
                  if ((drophel->Prob() > finehel->Prob()) && 
                      (ratdrop > 0.50) && (drophel->Fail() == 0)) {  //FIXME
                    if (4 == m_debug) cout << "append drop helix " <<  endl;
                    MdcxTrklist.append(drophel);
                    if(g_trklappend2) for (int iHit=0;iHit<drophel->Nhits();iHit++){g_trklappend2->fill(drophel->Layer(iHit));}
                    delete finehel;
                  } else {
                    if (4 == m_debug) cout << "append fine helix " <<  endl;
                    MdcxTrklist.append(finehel);
                    if(g_trklappend3)for(int iHit=0;iHit<finehel->Nhits();iHit++){g_trklappend3->fill(finehel->Layer(iHit));}
                    delete drophel;
                  }
                }

                //End hit dropping section
                int nwl = MdcxTrklist.length() - 1;
                if ((4 == m_debug) && (nwl > -1)) {
                  if (4 == m_debug) {
                    cout << endl << "found track " << trkCount<<std::endl;
                    MdcxTrklist[nwl]->print();
                    MdcxTrklist[nwl]->FitPrint();
                  }
                }
              } else { //*endif good fit so add to MdcxTrklist
                for (int i = 0; i < nSeg; i++) {
                  if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0);
                }
                delete finehel;
                trkCount--;
              }
            }//*end if not duplicate 
          }//*end floop first_prob >= probmin
        }// end iStV
      }// end iStU  
    }// end iAx
    if(4==m_debug) cout<< "end of this combo"<<endl;
  }// end iSegCombo
  if(4== m_debug) cout << " In MdcxFindTracks, found " << trkCount << " good tracks" << endl;
}//end of process

void MdcxFindTracks::resout

(

MdcxFittedHel *

f

)

[protected]

Definition at line 959 of file MdcxFindTracks.cxx.

References MdcxHel::Doca(), MdcxHel::Doca_Tof(), MdcxHel::Doca_Zh(), MdcxFittedHel::FitPrint(), m_debug, mhel, and MdcxFittedHel::XHitList().

                                                  {
  MdcxHel mdcxHel = (MdcxHel) *finehel;
  HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
  int kkk = 0; while(hitlist[kkk]){hitlist[kkk]->SetConstErr(0); kkk++;}
  MdcxFittedHel newhel(hitlist, mdcxHel);
  if(4 == m_debug) newhel.FitPrint(mhel,cout);
  kkk = 0;
  while(hitlist[kkk]) {
    float doca = newhel.Doca(*hitlist[kkk]);
    float zh   = newhel.Doca_Zh();
    float tof  = newhel.Doca_Tof();
    float dd   = hitlist[kkk]->d(newhel);
    float tt   = hitlist[kkk]->tcor(zh, tof);  
    if (4 == m_debug){
      cout << "MdcxFindTracks::resout (" 
        << hitlist[kkk]->Layer() << "," 
        << hitlist[kkk]->WireNo() << ") dt " 
        << tt << " resid " 
        << fabs(dd)-fabs(doca) << " doca " 
        << fabs(doca) << endl;
    }
    kkk++;
  }
  kkk=0; 
  while(hitlist[kkk]){hitlist[kkk]->SetConstErr(1); kkk++;}
}//endof drophits

MdcxHel MdcxFindTracks::TakeToOrigin

(

MdcxHel &

)

[protected]

Definition at line 986 of file MdcxFindTracks.cxx.

References MdcxHel::Code(), cos(), MdcxHel::CosPhi0(), MdcxHel::D0(), dlen(), MdcxHel::Omega(), MdcxHel::Phi0(), Constants::pi, sin(), MdcxHel::SinPhi0(), MdcxHel::T0(), MdcxHel::Tanl(), Constants::twoPi, MdcxHel::Xref(), MdcxHel::Yref(), and MdcxHel::Z0().

Referenced by process().

                                                  {
  double lng;
  double omegag = ihel.Omega();
  double phi0g = ihel.Phi0();
  double phi0s = phi0g; 
  double xp = ihel.Xref() - ihel.SinPhi0()*ihel.D0();
  double yp = ihel.Yref() + ihel.CosPhi0()*ihel.D0();
  double d0g = yp*cos(phi0g) - xp*sin(phi0g);
  if (omegag != 0.0) {
    phi0g = phi0g - Constants::pi;
    double xc = sin(phi0g)/omegag;
    double yc = -cos(phi0g)/omegag;
    double t1 = -xc - xp;
    double t2 = yc + yp;
    double phinew = atan2(t1, t2);
    if (omegag > 0.0) phinew += Constants::pi;
    if (phinew > Constants::pi) phinew -= Constants::twoPi;
    double xh = xc - sin(phinew)/omegag;
    double yh = yc + cos(phinew)/omegag;
    double x0 = xh + xp;
    double y0 = yh + yp;
    d0g = sqrt(x0*x0 + y0*y0);
    phi0g = phinew + Constants::pi;  
    double sd0 = x0*sin(phi0g) - y0*cos(phi0g);
    if (sd0 > 0.0) d0g = -d0g;
    lng = dlen(-2, phi0g, -1, phi0s, omegag);//yzhang TEMP 2011-10-17  
    //lng = dlen(phi0g, phi0s, omegag); 
  } else {
    lng = ihel.Xref()*ihel.CosPhi0() + ihel.Yref()*ihel.SinPhi0(); 
  }
  double tanlg = ihel.Tanl();
  double z0g = ihel.Z0() - lng*tanlg; 
  //cout << " z0g, tanlg, lng " << z0g << " " << tanlg << " " << lng << endl;
  double t0g = ihel.T0();
  int codeg = ihel.Code(); 
  MdcxHel ohel(d0g, phi0g, omegag, z0g, tanlg, t0g, codeg);
  return ohel;
}

bool MdcxFindTracks::testFromSameTrack	(	MdcxSeg *	seg1,
		MdcxSeg *	seg2
	)			[protected]

Definition at line 1025 of file MdcxFindTracks.cxx.

References genRecEmupikp::i, MdcxFittedHel::Nhits(), and MdcxFittedHel::XHitList().

Referenced by process().

                                                                  {
  /*
     std::cout<< " tkid seg1 "<<std::endl;
     for (int i =0; i<seg1->Nhits(); i++){
     int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
     std::cout<< tkid <<std::endl;
     }
     std::cout<< " tkid seg2 "<<std::endl;
     for (int i =0; i<seg2->Nhits(); i++){
     int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
     std::cout<< tkid <<std::endl;
     }
     */

  int trackId = -9999;
  for (int i =0; i<seg1->Nhits(); i++){
    int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
    if(tkid>=1000) tkid -=1000;
    if((trackId == -9999)){
      if(tkid>=0) trackId = tkid;
    }else{
      if(tkid!=trackId) return false;
    }
  } 
  for (int i =0; i<seg2->Nhits(); i++){
    int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
    if(tkid>=1000) tkid -=1000;
    if((trackId == -9999)){
      return false;
    }else{             
      if(tkid!=trackId) return false;
    }
  }
  return true;
}

---

Member Data Documentation

float MdcxFindTracks::curvmax [protected]

Definition at line 50 of file MdcxFindTracks.h.

Referenced by MdcxFindTracks().

int MdcxFindTracks::m_debug [protected]

Definition at line 57 of file MdcxFindTracks.h.

Referenced by drophits(), MdcxFindTracks(), process(), and resout().

bool MdcxFindTracks::m_doSag [protected]

Definition at line 58 of file MdcxFindTracks.h.

HepAList<MdcxFittedHel> MdcxFindTracks::MdcxTrklist [protected]

Definition at line 47 of file MdcxFindTracks.h.

Referenced by GetMdcxTrklist(), KillList(), and process().

MdcxHel MdcxFindTracks::mhel [protected]

Definition at line 48 of file MdcxFindTracks.h.

Referenced by process(), and resout().

int MdcxFindTracks::nbad [protected]

Definition at line 53 of file MdcxFindTracks.h.

Referenced by beginmore().

int MdcxFindTracks::nfail [protected]

Definition at line 54 of file MdcxFindTracks.h.

Referenced by beginmore().

int MdcxFindTracks::ngood [protected]

Definition at line 51 of file MdcxFindTracks.h.

Referenced by beginmore().

int MdcxFindTracks::nhitsgd [protected]

Definition at line 56 of file MdcxFindTracks.h.

Referenced by beginmore().

int MdcxFindTracks::nhitsmc [protected]

Definition at line 55 of file MdcxFindTracks.h.

Referenced by beginmore().

int MdcxFindTracks::npure [protected]

Definition at line 52 of file MdcxFindTracks.h.

Referenced by beginmore().

float MdcxFindTracks::probmin [protected]

Definition at line 49 of file MdcxFindTracks.h.

Referenced by MdcxFindTracks(), and process().

---