BOSS 8.0.0
BESIII Offline Software System
Loading...
Searching...
No Matches
MdcxFindTracks Class Reference

#include <MdcxFindTracks.h>

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
 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
 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)
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)
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< MdcxFittedHelMdcxTrklist
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

Constructor & Destructor Documentation

◆ MdcxFindTracks() [1/6]

◆ MdcxFindTracks() [2/6]

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

Definition at line 115 of file MdcxFindTracks.cxx.

115 {
117 curvmax = 1000000.0;
118 m_debug = debug;
119 process( segList );
120}
void process(const HepAList< MdcxSeg > &f)

◆ ~MdcxFindTracks() [1/3]

MdcxFindTracks::~MdcxFindTracks ( )
virtual

◆ MdcxFindTracks() [3/6]

MdcxFindTracks::MdcxFindTracks ( )

◆ MdcxFindTracks() [4/6]

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

◆ ~MdcxFindTracks() [2/3]

virtual MdcxFindTracks::~MdcxFindTracks ( )
virtual

◆ MdcxFindTracks() [5/6]

MdcxFindTracks::MdcxFindTracks ( )

◆ MdcxFindTracks() [6/6]

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

◆ ~MdcxFindTracks() [3/3]

virtual MdcxFindTracks::~MdcxFindTracks ( )
virtual

Member Function Documentation

◆ beginmore() [1/3]

◆ beginmore() [2/3]

void MdcxFindTracks::beginmore ( )
protected

◆ beginmore() [3/3]

void MdcxFindTracks::beginmore ( )
protected

◆ dlen() [1/3]

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

Definition at line 1053 of file MdcxFindTracks.cxx.

1053 {
1054 // double dp = p2 - p1;
1055 double dp; // yzhang change 2011-10-17
1056 if ( ( slay1 == 2 || slay1 == 3 || slay1 == 4 || slay1 == 9 || slay1 == 10 ) )
1057 { // Ax
1058 dp = p2 - p1;
1059 }
1060 else if ( ( slay2 >= slay1 ) ) { dp = p2 - p1; }
1061 else { dp = p1 - p2; }
1062
1063 if ( om < 0.0 )
1064 {
1065 while ( dp < -Constants::pi ) dp += Constants::pi;
1066 while ( dp > 0.0 ) dp -= 1 * Constants::pi;
1067 }
1068 else
1069 {
1070 while ( dp < 0.0 ) dp += 1 * Constants::pi;
1071 while ( dp > Constants::pi ) dp -= 1 * Constants::pi;
1072 }
1073 // double dl = dp * r;
1074 double dl = dp / om;
1075 // if(4 == m_debug) {
1076 // double r = 1./om;
1077 // cout prt(p1) prt(p2) prt(dp) prt(om) prt(r) prt(dl) << endl;
1078 // }
1079 return dl;
1080}
double p2[4]
double p1[4]

Referenced by process(), and TakeToOrigin().

◆ dlen() [2/3]

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

◆ dlen() [3/3]

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

◆ drophits() [1/3]

MdcxFittedHel MdcxFindTracks::drophits ( MdcxFittedHel * f)
protected

Definition at line 1082 of file MdcxFindTracks.cxx.

1082 {
1083 MdcxHel mdcxHel = (MdcxHel)*finehel;
1084 const HepAList<MdcxHit>& hitlist = finehel->XHitList();
1085 HepAList<MdcxHit> nwhitlist = hitlist;
1086 int nremove = 0;
1087 int kkk = 0;
1088 while ( nwhitlist[kkk] )
1089 {
1090 double pull = nwhitlist[kkk]->pull( mdcxHel );
1091 int layer = nwhitlist[kkk]->Layer();
1092
1093 if ( m_xtupleDropHits )
1094 {
1095 float t_doca = mdcxHel.Doca( *( nwhitlist[kkk] ) );
1096 float t_e = nwhitlist[kkk]->e( t_doca );
1098 m_segDropHitsLayer = layer;
1099 m_segDropHitsWire = nwhitlist[kkk]->WireNo();
1100 m_segDropHitsPull = pull;
1101 m_segDropHitsDoca = t_doca;
1102 m_segDropHitsSigma = t_e;
1103 m_segDropHitsDrift = nwhitlist[kkk]->d( mdcxHel );
1104 m_segDropHitsMcTkId = nwhitlist[kkk]->getDigi()->getTrackIndex();
1105 m_xtupleDropHits->write();
1106 }
1107
1108 if ( fabs( pull ) > MdcxParameters::dropHitsSigma[layer] )
1109 {
1110 nwhitlist.remove( kkk );
1111 nremove++;
1112 }
1113 else { kkk++; }
1114 }
1115 if ( m_debug == 4 )
1116 cout << " MdcxFindTracks::drophits drop " << nremove << " hits "
1117 << " nhits of drop helix = " << nwhitlist.length() << endl;
1118
1119 MdcxFittedHel newhel( nwhitlist, mdcxHel );
1120 return newhel;
1121} // endof drophits
int g_eventNo
Definition FTFinder.cxx:61
NTuple::Item< long > m_segDropHitsWire
NTuple::Item< double > m_segDropHitsDoca
NTuple::Item< double > m_segDropHitsSigma
NTuple::Item< long > m_segDropHitsLayer
NTuple::Item< double > m_segDropHitsMcTkId
NTuple::Item< long > m_segDropHitsEvtNo
NTuple::Item< double > m_segDropHitsPull
NTuple::Item< double > m_segDropHitsDrift
NTuple::Tuple * m_xtupleDropHits
double Doca(double WX, double WY, double WZ, double X, double Y, double Z=0.0)
Definition MdcxHel.cxx:261

Referenced by process().

◆ drophits() [2/3]

MdcxFittedHel MdcxFindTracks::drophits ( MdcxFittedHel * f)
protected

◆ drophits() [3/3]

MdcxFittedHel MdcxFindTracks::drophits ( MdcxFittedHel * f)
protected

◆ endmore() [1/3]

void MdcxFindTracks::endmore ( )
protected

Definition at line 1004 of file MdcxFindTracks.cxx.

1004{ print( cout ); } // endof endmore
void print(std::ostream &o)

◆ endmore() [2/3]

void MdcxFindTracks::endmore ( )
protected

◆ endmore() [3/3]

void MdcxFindTracks::endmore ( )
protected

◆ findz_cyl() [1/3]

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

Definition at line 1027 of file MdcxFindTracks.cxx.

1027 {
1028 double zint = -1000.0;
1029 double r = 1.0 / fabs( segList[iAx]->Omega() );
1030 double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0() * segList[iStU]->D0();
1031 double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0() * segList[iStU]->D0();
1032 double sx = segList[iStU]->Xline_slope();
1033 double sy = segList[iStU]->Yline_slope();
1034 double cx = xl_0 - segList[iAx]->Xc();
1035 double cy = yl_0 - segList[iAx]->Yc();
1036 double a = sx * sx + sy * sy;
1037 double b = 2.0 * ( sx * cx + sy * cy );
1038 // double c = (cx*cx + cy*cy - r*r);
1039 double bsq = b * b;
1040 double fac = 4.0 * a * ( cx * cx + cy * cy - r * r );
1041 double ta = 2.0 * a;
1042 if ( fac < bsq )
1043 {
1044 double sfac = sqrt( bsq - fac );
1045 double z1 = -( b - sfac ) / ta;
1046 double z2 = -( b + sfac ) / ta;
1047 zint = ( fabs( z1 ) < fabs( z2 ) ) ? z1 : z2;
1048 }
1049 // cout << " findz_cyl segs " << iAx << " " << iStU << " zint " << zint << endl;
1050 return zint;
1051}

Referenced by process().

◆ findz_cyl() [2/3]

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

◆ findz_cyl() [3/3]

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

◆ findz_sl() [1/3]

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

Definition at line 1006 of file MdcxFindTracks.cxx.

1006 {
1007 double Ap = -sin( segList[iAx]->Phi0_sl_approx() );
1008 double Bp = cos( segList[iAx]->Phi0_sl_approx() );
1009 double xp = segList[iAx]->Xref() + Ap * segList[iAx]->D0_sl_approx();
1010 double yp = segList[iAx]->Yref() + Bp * segList[iAx]->D0_sl_approx();
1011 double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0() * segList[iStU]->D0();
1012 double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0() * segList[iStU]->D0();
1013
1014 const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
1015 double Cl = hitsStU[0]->wz();
1016 double ndotl = Bp * hitsStU[0]->wy() + Ap * hitsStU[0]->wx();
1017 /*
1018 std::cout<<"ndotl "<<ndotl<<" "<<Bp<<" "<<hitsStU[0]->wy()<<" "<<
1019 Ap<<" "<<hitsStU[0]->wx()<<std::endl;
1020 */
1021 double zint =
1022 ( ndotl != 0.0 ) ? ( ( Bp * ( yp - yl ) + Ap * ( xp - xl ) ) * Cl / ndotl ) : -1000.;
1023 // cout << " findz_sl segs " << iAx << " " << iStU << " zint " << zint << endl;
1024 return zint;
1025}

Referenced by process().

◆ findz_sl() [2/3]

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

◆ findz_sl() [3/3]

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

◆ GetMdcxTrklist() [1/3]

const HepAList< MdcxFittedHel > & MdcxFindTracks::GetMdcxTrklist ( )
inline

◆ GetMdcxTrklist() [2/3]

const HepAList< MdcxFittedHel > & MdcxFindTracks::GetMdcxTrklist ( )
inline

◆ GetMdcxTrklist() [3/3]

const HepAList< MdcxFittedHel > & MdcxFindTracks::GetMdcxTrklist ( )
inline

◆ KillList() [1/3]

void MdcxFindTracks::KillList ( )
inlineprotected

Definition at line 60 of file InstallArea/x86_64-el9-gcc13-dbg/include/MdcxReco/MdcxFindTracks.h.

60{ HepAListDeleteAll( MdcxTrklist ); }

Referenced by ~MdcxFindTracks().

◆ KillList() [2/3]

void MdcxFindTracks::KillList ( )
inlineprotected

Definition at line 60 of file InstallArea/x86_64-el9-gcc13-opt/include/MdcxReco/MdcxFindTracks.h.

60{ HepAListDeleteAll( MdcxTrklist ); }

◆ KillList() [3/3]

void MdcxFindTracks::KillList ( )
inlineprotected

Definition at line 60 of file Reconstruction/MdcPatRec/MdcxReco/include/MdcxReco/MdcxFindTracks.h.

60{ HepAListDeleteAll( MdcxTrklist ); }

◆ plot() [1/3]

void MdcxFindTracks::plot ( ) const

Definition at line 986 of file MdcxFindTracks.cxx.

986{} // endof plot

◆ plot() [2/3]

void MdcxFindTracks::plot ( ) const

◆ plot() [3/3]

void MdcxFindTracks::plot ( ) const

◆ print() [1/3]

void MdcxFindTracks::print ( std::ostream & o)

Referenced by endmore().

◆ print() [2/3]

void MdcxFindTracks::print ( std::ostream & o)

◆ print() [3/3]

void MdcxFindTracks::print ( std::ostream & o)

◆ printpar() [1/3]

void MdcxFindTracks::printpar ( std::ostream & o)
protected

Referenced by beginmore().

◆ printpar() [2/3]

void MdcxFindTracks::printpar ( std::ostream & o)
protected

◆ printpar() [3/3]

void MdcxFindTracks::printpar ( std::ostream & o)
protected

◆ process() [1/3]

void MdcxFindTracks::process ( const HepAList< MdcxSeg > & f)

ZOUJH

ZOUJH

ZOUJH

ZOUJH

ZOUJH

Definition at line 124 of file MdcxFindTracks.cxx.

124 {
125 static int nfound = 0;
126 static float sumprob = 0.0;
127
128 MdcxHel thel( 0., 0., 0., 0., 0. );
129 mhel = thel;
130 int nSeg = segList.length();
131 for ( int i = 0; i < nSeg; i++ ) { segList[i]->SetUsedOnHel( 0 ); }
132
133 double xp, yp, xl1, yl1, rl1, xl2, yl2, rl2;
134 double d0g, phi0g, omegag, z0g, tanlg;
135 double d0g_sl, phi0g_sl, omegag_sl, z0g_sl, tanlg_sl;
136 double zintAU, zintAV, zintAU_sl, zintAV_sl = 9999.;
137 double rl1_sl, rl2_sl;
138 double xc = 1.e9, yc = 1.e9;// add initialize 25-05-15
139 double sinPhi0g_sl, cosPhi0g_sl, xp_sl, yp_sl;
140 double phiAx = 1.e9, phiStU, phiStV, xl1_sl, yl1_sl, xl2_sl, yl2_sl, xl1_0, yl1_0, xl2_0, yl2_0;// add initialize 25-05-15
141 double sx1, sy1, sx2, sy2;
142 double x0g, y0g;
143 double fltAx = 1.e9, ztest, fltL_sl, ztest_sl;// add initialize 25-05-15
144 int floop, trkCount = 0;
145
146 // The main loop; first get the 3 SL combo we'll use to make a trial helix
147 for ( int iSegCombo = 0; iSegCombo < MdcxParameters::nSegCombo; iSegCombo++ )
148 {
149 // FIXME optimize: group segList by superlayer ahead
150 int axlay = MdcxParameters::findTrkGroup[iSegCombo][0];
151 int stUlay = MdcxParameters::findTrkGroup[iSegCombo][1];
152 int stVlay = MdcxParameters::findTrkGroup[iSegCombo][2];
153
154 if ( 4 == m_debug )
155 std::cout << " Test combo: (" << axlay << "," << stUlay << "," << stVlay << ")"
156 << std::endl;
157
158 float max_dPhiAU = MdcxParameters::maxDp12[iSegCombo];
159 float max_dPhiAV = MdcxParameters::maxDp13[iSegCombo];
160
161 //---------------------------------------------------------------------
162 // Axial SL segment furnishes initial guess for d0, phi0, omega
163 //---------------------------------------------------------------------
164 for ( int iAx = 0; iAx < nSeg; iAx++ )
165 {
166 bool b_useGood_stU = true; // yzhang add 2009-10-21 17:20:45
167 bool b_useGood_stV = true; // yzhang add 2009-10-21 17:20:45
168 if ( ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iAx]->Fail() != 0 ) ||
169 ( segList[iAx]->SuperLayer( 0 ) != axlay ) || ( segList[iAx]->Origin() != -1 ) )
170 continue;
171
172 if ( 4 == m_debug )
173 {
174 std::cout << "---------1.Ax seg------ No." << iAx << std::endl;
175 segList[iAx]->printSegAll();
176 }
177
178 // Skip low pt track
179 omegag = segList[iAx]->Omega();
180 if ( fabs( omegag ) > MdcxParameters::maxTrkOmega )
181 { // FIXME 0.2
182 if ( 4 == m_debug ) std::cout << "SKIP by maxTrkOmega " << fabs( omegag ) << std::endl;
183 continue;
184 }
185 if ( 4 == m_debug ) std::cout << " iAx omegag = " << omegag << std::endl;
186 if ( g_omegag ) g_omegag->fill( omegag );
187
188 const HepAList<MdcxHit>& hitsAx = segList[iAx]->XHitList();
189
190 // Interpret segment as straight line also
191 d0g_sl = segList[iAx]->D0_sl_approx();
192 phi0g_sl = segList[iAx]->Phi0_sl_approx();
193 omegag_sl = 0.0;
194
195 // save c-tors; do it locally
196 sinPhi0g_sl = -sin( phi0g_sl );
197 cosPhi0g_sl = cos( phi0g_sl );
198 xp_sl = segList[iAx]->Xref() + sinPhi0g_sl * d0g_sl;
199 yp_sl = segList[iAx]->Yref() + cosPhi0g_sl * d0g_sl;
200 d0g_sl = yp_sl * cosPhi0g_sl + xp_sl * sinPhi0g_sl;
201 fltL_sl = segList[iAx]->Xref() * cosPhi0g_sl - segList[iAx]->Yref() * sinPhi0g_sl;
202 if ( 4 == m_debug )
203 {
204 std::cout << "--Ax :" prt( d0g_sl ) prt( phi0g_sl ) prt( omegag_sl ) prt( sinPhi0g_sl )
205 prt( xp_sl ) prt( yp_sl ) prt( fltL_sl )
206 << std::endl;
207 }
208
209 if ( fabs( omegag ) > Constants::epsilon )
210 {
211 MdcxHel ohel = TakeToOrigin( *segList[iAx] );
212 omegag = ohel.Omega();
213 phi0g = ohel.Phi0();
214 d0g = ohel.D0();
215 xc = ohel.Xc();
216 yc = ohel.Yc();
217 x0g = ohel.X0();
218 y0g = ohel.Y0();
219 phiAx = atan2( y0g - yc, x0g - xc );
220 fltAx = dlen( -2, phi0g, -1, segList[iAx]->Phi0(), omegag );
221 // fltAx = dlen(phi0g,segList[iAx]->Phi0(), omegag); //yzhang TEMP 2011-10-17
222 if ( 4 == m_debug )
223 {
224 std::cout << "--Ax :"
225 << " ohel ";
226 ohel.print();
227 }
228 }
229
230 //---------------------------------------------------------------------
231 // First stereo SL segement and axial segment furnish one z intersection
232 //---------------------------------------------------------------------
233 // for (int iStU = 0; iStU < nSeg; iStU++ ) //yzhang delete 2009-10-21 18:13:50
234 for ( int iStU = -1; true; )
235 { // yzhang add
236 if ( !b_useGood_stU && ( iStU >= ( nSeg - 1 ) ) ) break;
237 if ( b_useGood_stU && ( iStU >= ( nSeg - 1 ) ) )
238 {
239 iStU = -1;
240 b_useGood_stU = false;
241 }
242 iStU++;
243 if ( ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iStU]->GetUsedOnHel() != 0 ) ||
244 ( segList[iStU]->Fail() != 0 ) || ( segList[iStU]->SuperLayer( 0 ) != stUlay ) ||
245 ( segList[iStU]->Origin() != -1 )
246 // yzhang add 2009-10-21 17:21:55
247 || ( b_useGood_stU && ( segList[iStU]->Nhits() < 4 ) ) )
248 { continue; }
249
250 if ( 4 == m_debug )
251 {
252 std::cout << " Test combo: AUV " << axlay << "," << stUlay << "," << stVlay
253 << std::endl;
254 std::cout << "---------2.St U seg ------No." << iStU << std::endl; // yzhang debug
255 segList[iStU]->printSeg();
256 std::cout << " omega of seg iStU = " << segList[iStU]->Omega() << std::endl;
257 }
258
259 // Test dPhi between Ax slayer and StU slayer
260 const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
261 double dPhiAU = fabs( hitsAx[0]->phiMid() - hitsStU[0]->phiMid() );
262 if ( dPhiAU > Constants::pi ) dPhiAU = Constants::twoPi - dPhiAU;
263 if ( 4 == m_debug )
264 {
265 if ( dPhiAU > max_dPhiAU )
266 { cout << "SKIP by dPhiAU " << dPhiAU << " > " << max_dPhiAU << endl; }
267 else { cout << "KEEP by dPhiAU " << dPhiAU << " <= " << max_dPhiAU << endl; }
268 std::cout << " phi mid Ax " << hitsAx[0]->phiMid() << " StU " << hitsStU[0]->phiMid()
269 << std::endl;
270 std::cout << std::endl;
271 }
272 if ( g_dPhiAU ) g_dPhiAU->fill( dPhiAU, segList[iStU]->SuperLayer() );
273
274 if ( dPhiAU > max_dPhiAU ) { continue; }
275
276 // Calculate z by Ax and StU segs
277 xl1_0 = segList[iStU]->Xline_bbrrf();
278 yl1_0 = segList[iStU]->Yline_bbrrf();
279 sx1 = segList[iStU]->Xline_slope();
280 sy1 = segList[iStU]->Yline_slope();
281 if ( 4 == m_debug )
282 std::cout prt( xl1_0 ) prt( yl1_0 ) prt( sx1 ) prt( sy1 ) prt( omegag ) << std::endl;
283 if ( fabs( omegag ) > Constants::epsilon )
284 {
285 zintAU = findz_cyl( iAx, iStU, segList );
286 xl1 = xl1_0 + sx1 * zintAU;
287 yl1 = yl1_0 + sy1 * zintAU;
288 rl1 = sqrt( xl1 * xl1 + yl1 * yl1 );
289 phiStU = atan2( yl1 - yc, xl1 - xc );
290 if ( 4 == m_debug )
291 cout prt( zintAU ) prt( xl1 ) prt( yl1 ) prt( rl1 ) prt( phiStU ) << endl;
292 }
293 else { zintAU = -9999.; }
294
295 // Calculate z_sl by Ax and StU segs
296 zintAU_sl = findz_sl( iAx, iStU, segList );
297 xl1_sl = xl1_0 + sx1 * zintAU_sl;
298 yl1_sl = yl1_0 + sy1 * zintAU_sl;
299 rl1_sl = sqrt( xl1_sl * xl1_sl + yl1_sl * yl1_sl );
300 if ( 4 == m_debug )
301 cout prt( zintAU_sl ) prt( xl1_sl ) prt( yl1_sl ) prt( rl1_sl ) << endl;
302
303 // Cut by z and z_sl
304 if ( ( zintAU < -MdcxParameters::maxMdcZLen ) &&
305 ( zintAU_sl < -MdcxParameters::maxMdcZLen ) )
306 {
307 if ( 4 == m_debug )
308 std::cout << " SKIP by zintAU < max MdcZLen " << MdcxParameters::maxMdcZLen
309 << std::endl;
310 continue;
311 }
312 if ( ( zintAU > MdcxParameters::maxMdcZLen ) &&
313 ( zintAU_sl > MdcxParameters::maxMdcZLen ) )
314 {
315 if ( 4 == m_debug )
316 std::cout << " SKIP by zintAU > max MdcZLen " << MdcxParameters::maxMdcZLen
317 << std::endl;
318 continue;
319 }
320
321 //---------------------------------------------------------------------
322 // Next stereo SL segement and axial segment furnish other z intersection
323 //---------------------------------------------------------------------
324 for ( int iStV = -1; true; )
325 {
326 if ( !b_useGood_stV && ( iStV >= ( nSeg - 1 ) ) )
327 {
328 if ( 4 == m_debug ) std::cout << iStV << " no segments in V slayer " << std::endl;
329 break;
330 }
331 if ( b_useGood_stV && ( iStV >= ( nSeg - 1 ) ) )
332 {
333 iStV = -1;
334 b_useGood_stV = false;
335 }
336 iStV++;
337 if ( ( segList[iStV]->GetUsedOnHel() != 0 ) ||
338 ( segList[iStU]->GetUsedOnHel() != 0 ) ||
339 ( segList[iAx]->GetUsedOnHel() != 0 ) || ( segList[iStV]->Fail() != 0 ) ||
340 ( segList[iStV]->SuperLayer( 0 ) != stVlay ) ||
341 ( segList[iStV]->Origin() != -1 )
342 // yzhang add 2009-10-21 17:21:55
343 || ( b_useGood_stV && ( segList[iStV]->Nhits() < 4 ) ) )
344 { continue; }
345
346 if ( 4 == m_debug )
347 {
348 std::cout << " Test combo: AUV " << axlay << "," << stUlay << "," << stVlay
349 << std::endl;
350 std::cout << "---------3.St V seg ------No." << iStV << std::endl; // yzhang debug
351 segList[iStV]->printSeg();
352 std::cout << " omega of seg iStV = " << segList[iStV]->Omega() << std::endl;
353 }
354
355 // Calculate dPhi between Ax and StV segs
356 const HepAList<MdcxHit>& hitsStV = segList[iStV]->XHitList();
357 double dPhiAV = fabs( hitsAx[0]->phiMid() - hitsStV[0]->phiMid() );
358 if ( dPhiAV > Constants::pi ) dPhiAV = Constants::twoPi - dPhiAV;
359
360 if ( 4 == m_debug )
361 {
362 if ( dPhiAV > max_dPhiAV )
363 { cout << "SKIP by dPhiAV " << dPhiAV << " > " << max_dPhiAV << endl; }
364 else { cout << "KEEP by dPhiAV " << dPhiAV << " <= " << max_dPhiAV << endl; }
365 std::cout << " phi mid Ax " << hitsAx[0]->phiMid() << " StV "
366 << hitsStV[0]->phiMid() << std::endl;
367 std::cout << std::endl;
368 }
369 if ( g_dPhiAV ) g_dPhiAV->fill( dPhiAV, segList[iStV]->SuperLayer() );
370 if ( dPhiAV > max_dPhiAV ) { continue; }
371
372 // Calculate z by Ax and StV segs
373 xl2_0 = segList[iStV]->Xline_bbrrf();
374 yl2_0 = segList[iStV]->Yline_bbrrf();
375 sx2 = segList[iStV]->Xline_slope();
376 sy2 = segList[iStV]->Yline_slope();
377 if ( fabs( omegag ) > Constants::epsilon )
378 {
379 zintAV = findz_cyl( iAx, iStV, segList );
380 xl2 = xl2_0 + sx2 * zintAV;
381 yl2 = yl2_0 + sy2 * zintAV;
382 rl2 = sqrt( xl2 * xl2 + yl2 * yl2 );
383 if ( 4 == m_debug )
384 {
385 segList[iAx]->printSeg();
386 segList[iStV]->printSeg();
387 cout << "zintAV " << zintAV << endl;
388 }
389 phiStV = atan2( yl2 - yc, xl2 - xc );
390 }
391 else { zintAV = -9999.; }
392
393 // Calculate z by Ax and StV segs
394 zintAV_sl = findz_sl( iAx, iStV, segList );
395 xl2_sl = xl2_0 + sx2 * zintAV_sl;
396 yl2_sl = yl2_0 + sy2 * zintAV_sl;
397 rl2_sl = sqrt( xl2_sl * xl2_sl + yl2_sl * yl2_sl );
398
399 // Cut by z of AV segs
400 if ( ( zintAV < -MdcxParameters::maxMdcZLen ) &&
401 ( zintAV_sl < -MdcxParameters::maxMdcZLen ) )
402 {
403 if ( 4 == m_debug )
404 std::cout << " SKIP by zintAV < max MdcZLen " << MdcxParameters::maxMdcZLen
405 << std::endl;
406 continue;
407 }
408 if ( ( zintAV > MdcxParameters::maxMdcZLen ) &&
409 ( zintAV_sl > MdcxParameters::maxMdcZLen ) )
410 {
411 if ( 4 == m_debug )
412 std::cout << " SKIP by zintAV > max MdcZLen " << MdcxParameters::maxMdcZLen
413 << std::endl;
414 continue;
415 }
416
417 // We now have enough info to form and fit a guess helix (inside iAx-iStU-iStV loops)
418 if ( 4 == m_debug )
419 {
420 // cout << "KEEP layer[klm]: [" << axlay <<","<< stUlay ","<< stVlay << "] "
421 // prt(d0g) prt(phi0g) prt(omegag) zrt(zintAU) prt(zintAV) << endl;
422 segList[iAx]->printSeg();
423 segList[iStU]->printSeg();
424 segList[iStV]->printSeg();
425 }
426 MdcxFittedHel fithel;
427 double helixFitSigma = MdcxParameters::helixFitSigma;
428 double first_prob = 0.0;
429 HepAList<MdcxHit> hitlist;
430 hitlist.append( hitsAx );
431 hitlist.append( hitsStU );
432 hitlist.append( hitsStV );
433 if ( g_trkllmk )
434 for ( int i = 0; i < hitlist.length(); i++ )
435 { g_trkllmk->fill( hitlist[i]->Layer() ); }
436
437 //---------------------------------------------------------------------
438 // Try fitting the trial helix from axial SL as circle first
439 //---------------------------------------------------------------------
440 floop = 1;
441 while ( floop )
442 {
443 if ( 4 == m_debug )
444 std::cout << "---------4.Ax circle fit------" << std::endl; // yzhang debug
445 if ( fabs( omegag ) < Constants::epsilon )
446 {
447 if ( 4 == m_debug ) std::cout << "SKIP by omegag==0 " << std::endl;
448 break;
449 }
450 if ( ( zintAU < -MdcxParameters::maxMdcZLen ) ||
451 ( zintAU > MdcxParameters::maxMdcZLen ) )
452 {
453 if ( 4 == m_debug )
454 std::cout << "SKIP by zintAU out of range " << zintAU << " "
455 << MdcxParameters::maxMdcZLen << std::endl;
456 break;
457 }
458 if ( ( zintAV < -MdcxParameters::maxMdcZLen ) ||
459 ( zintAV > MdcxParameters::maxMdcZLen ) )
460 {
461 if ( 4 == m_debug )
462 std::cout << "SKIP by zintAU out of range " << zintAU << " "
463 << MdcxParameters::maxMdcZLen << std::endl;
464 break;
465 }
466
467 // Calculate flight length of AUV segs
468 if ( 4 == m_debug )
469 cout << "dlen calc. slay U " << segList[iStU]->SuperLayer() << " slay V "
470 << segList[iStV]->SuperLayer() << endl;
471 double fltLenUV =
472 dlen( segList[iStU]->SuperLayer(), phiStU, segList[iStV]->SuperLayer(), phiStV,
473 omegag ); // yzhang change TEMP 2011-10-17
474 // double fltLenUV = dlen(phiStU, phiStV, omegag);
475 if ( 4 == m_debug )
476 { cout prt( fltLenUV ) prt( phiStU ) prt( phiStV ) prt( omegag ) << std::endl; }
477 if ( fltLenUV > MdcxParameters::maxDlen )
478 { // 15. to 150.
479 if ( 4 == m_debug )
480 std::cout << "SKIP by fltLenUV > " << MdcxParameters::maxDlen << std::endl;
481 break; // FIXME
482 }
483 tanlg = ( zintAV - zintAU ) / fltLenUV;
484 if ( 4 == m_debug )
485 cout << "dlen calc. slay A " << segList[iAx]->SuperLayer() << " slay U "
486 << segList[iStU]->SuperLayer() << endl;
487 double fltLenAU =
488 dlen( segList[iAx]->SuperLayer(), phiAx, segList[iStU]->SuperLayer(), phiStU,
489 omegag ); // yzhang change TEMP 2011-10-17
490
491 if ( m_xtupleSegComb )
492 {
493 m_segCombOmega = omegag;
494 m_segCombSameAU = testFromSameTrack( segList[iAx], segList[iStU] );
495 m_segCombSameUV = testFromSameTrack( segList[iStU], segList[iStV] );
496 m_segCombSlA = segList[iAx]->SuperLayer();
497 m_segCombSlU = segList[iStU]->SuperLayer();
498 m_segCombSlV = segList[iStV]->SuperLayer();
499 m_segCombPhiA = phiAx;
500 m_segCombPhiU = phiStU;
501 m_segCombPhiV = phiStV;
502 m_segCombDLenUV = fltLenUV;
503 m_segCombDLenAU = fltLenAU;
504 m_xtupleSegComb->write();
505 }
506 // double fltLenAU = dlen(phiAx, phiStU, omegag);
507 z0g = zintAU - fltLenAU * tanlg;
508 if ( 4 == m_debug )
509 cout prt( z0g ) prt( tanlg ) prt( fltLenAU ) prt( zintAU ) prt( zintAV ) << endl;
510
511 // Cut by z0g
512 if ( ( z0g < -MdcxParameters::maxMdcZLen ) ||
514 {
515 if ( 4 == m_debug )
516 std::cout << "SKIP by z0g out of range " << z0g << ">"
517 << MdcxParameters::maxMdcZLen << std::endl;
518 break;
519 }
520 ztest = z0g + fltAx * tanlg;
521 if ( 4 == m_debug ) std::cout prt( ztest ) prt( fltAx ) << std::endl;
522
523 //---Helix fitting
524 // MdcxHel ghel(d0g,phi0g,omegag,z0g,tanlg); // ghel.print();
525 MdcxHel ghel( segList[iAx]->D0(), segList[iAx]->Phi0(), segList[iAx]->Omega(),
526 ztest, tanlg, 0.0, 11111, 0, segList[iAx]->Xref(),
527 segList[iAx]->Yref() );
528 // zoujh?: ghel.SetTurnFlag(1); // ghel.print();
529 MdcxFittedHel firstfit( hitlist, ghel, helixFitSigma );
530 first_prob = firstfit.Prob();
531 if ( 4 == m_debug )
532 {
533 std::cout << " after first fit: " << std::endl;
534 firstfit.FitPrint();
535 }
536 if ( first_prob >= probmin )
537 {
538 MdcxHel helixOrigin = TakeToOrigin( firstfit );
539 MdcxFittedHel nextfit( hitlist, helixOrigin, helixFitSigma );
540 first_prob = nextfit.Prob();
541 fithel = nextfit;
542 if ( g_trklgood )
543 for ( int i = 0; i < nextfit.Nhits(); i++ )
544 { g_trklgood->fill( nextfit.Layer( i ) ); }
545 }
546 if ( g_trklfirstProb ) g_trklfirstProb->fill( first_prob );
547 if ( 4 == m_debug ) cout << " prob after helix fitting = " << first_prob << endl;
548 floop = 0;
549 } // end fitting trial helix
550
551 //---------------------------------------------------------------------
552 // Try fitting straight line trial helix if fit from axial-SL-as-circle trial
553 // if helix not good enough
554 //---------------------------------------------------------------------
555 floop = 1;
556 while ( floop )
557 {
558 if ( 4 == m_debug )
559 std::cout << "---------4.2 straight line fit------" << std::endl; // yzhang debug
560 if ( 4 == m_debug )
561 cout << " zintAU_sl " << zintAU_sl << " zintAV_sl " << zintAV_sl << endl;
562 if ( first_prob > probmin )
563 {
564 if ( 4 == m_debug ) cout << " helix fit good , exit straight line fit." << endl;
565 break;
566 }
567
568 if ( ( zintAU_sl < -MdcxParameters::maxMdcZLen ) ||
569 ( zintAU_sl > MdcxParameters::maxMdcZLen ) )
570 break;
571 if ( ( zintAV_sl < -MdcxParameters::maxMdcZLen ) ||
572 ( zintAV_sl > MdcxParameters::maxMdcZLen ) )
573 break;
574
575 double dx = xl2_sl - xl1_sl;
576 double dy = yl2_sl - yl1_sl;
577 double dl = sqrt( dx * dx + dy * dy );
578 tanlg_sl = ( zintAV_sl - zintAU_sl ) / dl;
579 dx = xl1_sl + d0g_sl * sin( phi0g_sl );
580 dy = yl1_sl - d0g_sl * cos( phi0g_sl );
581 dl = sqrt( dx * dx + dy * dy );
582 z0g_sl = zintAU_sl - dl * tanlg_sl;
583 if ( 4 == m_debug )
584 cout prt( d0g_sl ) prt( phi0g_sl ) prt( z0g_sl ) prt( tanlg_sl ) << endl;
585
586 if ( ( z0g_sl < -MdcxParameters::maxMdcZLen ) ||
587 ( z0g_sl > MdcxParameters::maxMdcZLen ) )
588 {
589 if ( 4 == m_debug )
590 std::cout << "SKIP by z0g_sl:" << z0g_sl << " > " << MdcxParameters::maxMdcZLen
591 << std::endl;
592 break;
593 }
594 // MdcxHel ghel(d0g_sl,phi0g_sl,omegag_sl,z0g_sl,tanlg_sl); // ghel.print();
595 ztest_sl = z0g_sl + fltL_sl * tanlg_sl;
596 if ( 4 == m_debug ) cout prt( ztest_sl ) << endl;
597
598 // Helix fitting
599 MdcxHel ghel( segList[iAx]->D0_sl_approx(), segList[iAx]->Phi0_sl_approx(),
600 omegag_sl, ztest_sl, tanlg_sl, 0.0, 11111, 0, segList[iAx]->Xref(),
601 segList[iAx]->Yref() );
602 // zoujh?: ghel.SetTurnFlag(1); // ghel.print();
603 MdcxFittedHel firstfit( hitlist, ghel, helixFitSigma );
604 first_prob = firstfit.Prob();
605 // if (first_prob >= probmin)fithel=firstfit;
606 if ( 4 == m_debug )
607 {
608 ghel.print();
609 std::cout << "first_prob " << first_prob << std::endl;
610 }
611
612 if ( first_prob >= probmin )
613 {
614 MdcxHel helixOrigin = TakeToOrigin( firstfit );
615 MdcxFittedHel nextfit( hitlist, helixOrigin, helixFitSigma );
616 first_prob = nextfit.Prob();
617 fithel = nextfit;
618 if ( 4 == m_debug )
619 {
620 cout << " three seg sl nextfit" << endl;
621 nextfit.FitPrint();
622 }
623 }
624 floop = 0;
625 }
626
627 //-----------------------------------------------------------------
628 // Got a good trial helix fit; now try to add other segments onto it
629 //-----------------------------------------------------------------
630 floop = 1;
631 while ( floop )
632 {
633 floop = 0;
634 if ( 4 == m_debug )
635 std::cout << "---------5. try add seg to helix------" << std::endl;
636 if ( first_prob < probmin )
637 {
638 if ( 4 == m_debug )
639 std::cout << "prob" << first_prob << " " << probmin
640 << " fit NOT good, exit add segs " << std::endl;
641 break;
642 }
643 float bchisq = fithel.Chisq() * helixFitSigma * helixFitSigma;
644 int bndof = fithel.Nhits() - 5;
645 float bprob = Mdcxprobab( bndof, bchisq );
646 trkCount++;
647 segList[iAx]->SetUsedOnHel( trkCount );
648 segList[iStU]->SetUsedOnHel( trkCount );
649 segList[iStV]->SetUsedOnHel( trkCount );
650
651 if ( 4 == m_debug )
652 {
653 cout << "in add seg to trail helix, klm seg:" << endl; /// ZOUJH
654 segList[iAx]->printSeg(); /// ZOUJH
655 segList[iStU]->printSeg(); /// ZOUJH
656 segList[iStV]->printSeg(); /// ZOUJH
657 }
658
659 // Loop superlayers
660 for ( int iSlay = 0; iSlay < Constants::nSuperLayer; iSlay++ )
661 {
662 int ilay = MdcxParameters::layerSet2AddSeg[iSegCombo][iSlay];
663 for ( int iSeg = 0; iSeg < nSeg; iSeg++ )
664 {
665 if ( ( segList[iSeg]->SuperLayer( 0 ) != ilay ) ||
666 ( segList[iSeg]->GetUsedOnHel() != 0 ) ||
667 ( segList[iSeg]->Fail() != 0 ) || ( segList[iSeg]->Origin() != -1 ) )
668 continue;
669 if ( 4 == m_debug )
670 {
671 std::cout << endl << "== add seg " << iSeg << " ";
672 segList[iSeg]->printSeg();
673 }
674
675 // Cut by phi_diff = (Phi_Ax - Phi_Add)
676 float phiAxSeg = segList[iAx]->XHitList()[0]->phiMid();
677 float phiAddSeg = segList[iSeg]->XHitList()[0]->phiMid();
678 double phiKNSegDiff = fabs( phiAxSeg - phiAddSeg );
679 // yzhang add 2011-10-11
680 // double phiKNSegDiff = (phiAxSeg - phiAddSeg);
681 // if(phiKNSegDiff<-Constants::pi) phiKNSegDiff+=2*Constants::pi;
682 // if(phiKNSegDiff>Constants::pi) phiKNSegDiff-=2*Constants::pi;
683 int t_same = -999;
684 if ( m_xtupleAddSeg1 )
685 {
686 // std::cout<< __FILE__ << " testFromSameTrack " << fromSameTrack<< "
687 // "<<std::endl;
688 t_same = testFromSameTrack( segList[iAx], segList[iSeg] );
689 m_addSegSame = t_same;
690 m_addSegSeedSl = segList[iAx]->SuperLayer();
691 m_addSegSeedPhi = segList[iAx]->XHitList()[0]->phiMid();
692 m_addSegSeedPhiLay = segList[iAx]->XHitList()[0]->Layer();
693 m_addSegSeedPhi0 = segList[iAx]->Phi0_sl_approx();
694 m_addSegSeedD0 = segList[iAx]->D0_sl_approx();
695 m_addSegAddSl = segList[iSeg]->SuperLayer();
696 m_addSegAddPhi = segList[iSeg]->XHitList()[0]->phiMid();
697 m_addSegAddPhiLay = segList[iSeg]->XHitList()[0]->Layer();
698 m_addSegAddPhi0 = segList[iSeg]->Phi0_sl_approx();
699 m_addSegAddD0 = segList[iSeg]->D0_sl_approx();
700 m_xtupleAddSeg1->write();
701 }
702 // yzhang change 2011-10-11
703 if ( phiKNSegDiff > 0.5 * Constants::pi && phiKNSegDiff < 1.5 * Constants::pi )
704 // if (fabs(phiKNSegDiff) > 0.7*Constants::pi )
705 { // yzhang TEMP
706 if ( 4 == m_debug )
707 std::cout << " SKIP by phi diff,same " << t_same << " Ax " << phiAxSeg
708 << " iSeg " << phiAddSeg << " diff " << phiKNSegDiff
709 << std::endl; // yzhang debug
710 continue; // zoujh FIXME
711 }
712 else
713 {
714 if ( 4 == m_debug )
715 std::cout << " phi diff OK, Ax " << phiAxSeg << " added " << phiAddSeg
716 << " diff=" << phiKNSegDiff << std::endl; // yzhang debug
717 }
718
719 // Calculate doca
720 xp = segList[iSeg]->Xref() - segList[iSeg]->SinPhi0() * segList[iSeg]->D0();
721 yp = segList[iSeg]->Yref() + segList[iSeg]->CosPhi0() * segList[iSeg]->D0();
722 const HepAList<MdcxHit>& hitsSegAdd = segList[iSeg]->XHitList();
723 double proca = fithel.Doca( hitsSegAdd[0]->wx(), hitsSegAdd[0]->wy(),
724 hitsSegAdd[0]->wz(), xp, yp );
725 if ( g_trklproca ) g_trklproca->fill( proca );
726
727 // Segment passes loose doca cut; does it belong on this track?
728 if ( m_xtupleAddSeg2 )
729 {
731 m_addSegPoca = proca;
732 m_addSegSlayer = iSlay;
733 m_addSegLen = fithel.Doca_Len();
734 m_xtupleAddSeg2->write();
735 }
736
737 if ( !( ( fabs( proca ) < MdcxParameters::maxProca ) &&
738 ( fithel.Doca_Len() < fithel.Lmax() ) ) )
739 {
740 // if (g_trklprocaSl) g_trklprocaSl->fill(segList[iSeg]->SuperLayer(0));
741 if ( 4 == m_debug )
742 {
743 std::cout << " SKIP by fabs(proca):" << fabs( proca ) << "<"
745 << " or Doca_Len:" << fithel.Doca_Len() << "<" << fithel.Lmax()
746 << std::endl; // yzhang debug
747 }
748 }
749 else
750 {
751 if ( 4 == m_debug )
752 {
753 std::cout << " proca & len OK, |proca| " << fabs( proca ) << "<"
755 << " && Doca_Len:" << fithel.Doca_Len() << "<" << fithel.Lmax()
756 << std::endl; // yzhang debug
757 }
758 }
759 if ( fithel.Doca_Len() < 0 ) continue; // yzhang add TEMP 2011-10-11
760
761 if ( ( fabs( proca ) < MdcxParameters::maxProca ) &&
762 ( fithel.Doca_Len() < fithel.Lmax() ) )
763 {
764 MdcxFittedHel newhel;
765 newhel.Grow( fithel, hitsSegAdd );
766 if ( newhel.Nhits() != hitlist.length() )
767 {
768 if ( 4 == m_debug )
769 {
770 cout << " rcs " << newhel.Rcs() << "<=?"
772 }
773 // yzhang change 2009-10-21 FIXME
774 // if (newhel.Rcs() > MdcxParameters::maxRcsInAddSeg)
775 if ( newhel.Prob() < probmin )
776 {
777 if ( 4 == m_debug )
778 {
779 cout << " prob " << newhel.Prob() << "<" << probmin << ", ReFit"
780 << endl; /// ZOUJH
781 }
782 newhel.ReFit(); // FIXME
783 }
784 if ( 4 == m_debug )
785 {
786 cout << " refit prob " << newhel.Prob() << "<" << probmin << " Fail? "
787 << newhel.Fail() << endl;
788 }
789 // yzhang change 2009-10-21 FIXME
790 if ( ( newhel.Prob() >= probmin ) && ( newhel.Fail() == 0 ) )
791 {
792 // if ( (newhel.Rcs() <= MdcxParameters::maxRcsInAddSeg) &&
793 // (newhel.Fail() == 0) )
794 fithel = newhel;
795 hitlist = newhel.XHitList();
796 segList[iSeg]->SetUsedOnHel( trkCount );
797 bchisq = fithel.Chisq() * helixFitSigma * helixFitSigma;
798 bndof = fithel.Nhits() - 5;
799 float tprob = Mdcxprobab( bndof, bchisq );
800 if ( tprob > bprob ) bprob = tprob;
801 if ( 4 == m_debug )
802 { cout << " segment ADDED, prob=" << newhel.Prob() << endl; }
803 } //*end trial track still good so keep new seg
804 else
805 {
806 if ( 4 == m_debug ) std::cout << " fit FAILED" << std::endl;
807 }
808 }
809 else
810 { //*end new hits so refit
811 segList[iSeg]->SetUsedOnHel( trkCount );
812 if ( 4 == m_debug ) cout << " segment ADDED, but no new hits" << endl;
813 } //*end dont need to refit if no new hits
814 } //*end fabs(proca) < 0.010
815 } //*end iSeg (loop over segs in a superlayer)
816 } //*end iSlay (loop over superlayers)
817 if ( ( ( fithel.Prob() < 0.90 ) && ( fithel.Nhits() == 12 ) ) ||
818 ( fithel.Fail() != 0 ) )
819 {
820 for ( int i = 0; i < nSeg; i++ )
821 {
822 if ( segList[i]->GetUsedOnHel() == trkCount ) segList[i]->SetUsedOnHel( 0 );
823 }
824 trkCount--;
825 break; //*jump to end floop first_prob >= probmin
826 }
827
828 if ( 4 == m_debug )
829 {
830 std::cout << "Track after add segs" << std::endl; // yzhang debug
831 fithel.FitPrint();
832 }
833
834 // See if helix better moving backwards
835 if ( 4 == m_debug )
836 std::cout << "---------6. flip------" << std::endl; // yzhang debug
837 fithel.flip();
838
839 // See if this track shares lots of hits with another saved track
840 int kki = 0;
841 int notduplicate = 1;
842 while ( MdcxTrklist[kki] )
843 {
844 if ( 4 == m_debug )
845 std::cout << "---------7. test saved track------" << std::endl; // yzhang debug
846 const HepAList<MdcxHit>& junk = MdcxTrklist[kki]->XHitList();
847 int overl = 0;
848 int kkj = 0;
849 while ( junk[kkj] )
850 {
851 int kkl = 0;
852 while ( hitlist[kkl] )
853 {
854 if ( hitlist[kkl] == junk[kkj] ) overl++;
855 kkl++;
856 }
857 kkj++;
858 }
859 if ( 4 == m_debug )
860 cout << "overlap with track # " << kki << " is " << junk.length()
861 << " hitList " << hitlist.length() << " overlap " << overl << endl;
862 if ( ( hitlist.length() / 4. ) <= overl ) notduplicate = 0;
863 // yzhang 2009-10-21 16:30:53
864 // if ( hitlist.length() <= (overl+2) ) notduplicate = 0;
865 kki++;
866 }
867
868 if ( g_trklhelix )
869 for ( int iHit = 0; iHit < fithel.Nhits(); iHit++ )
870 { g_trklhelix->fill( fithel.Layer( iHit ) ); }
871
872 // If not a duplicate, is it worth saving?
873 if ( notduplicate )
874 {
875 if ( 4 == m_debug )
876 std::cout << "---------8. test worth saving?------"
877 << std::endl; // yzhang debug
878 MdcxFittedHel* finehel = new MdcxFittedHel( hitlist, fithel );
879
880 if ( 4 == m_debug )
881 {
882 std::cout << __FILE__ << " finehel Prob>0.001 " << finehel->Prob() << " nhits "
883 << finehel->Nhits() << ">=15? "
884 << " bprob " << bprob << " >=?probmin " << probmin << " Fail==0? "
885 << finehel->Fail() << std::endl;
886 finehel->FitPrint();
887 }
888 // cut nhits from 15->9 yzhang 2009-10-21
889 if ( ( ( finehel->Prob() > 0.001 ) || ( finehel->Nhits() >= 15 ) ||
890 ( bprob > probmin )
891 //|| (bchisq/bndof > MdcxParameters::maxRcsInAddSeg)
892 ) &&
893 ( finehel->Fail() == 0 ) )
894 {
895 nfound++;
896 sumprob += finehel->Prob();
897 // if ( hitlist.length() == 16 )resout(finehel);
898
899 // Begin hit dropping section
900 // Set nodrop != 0 if want to prohibit dropping bad hits
901 int nodrop = 0;
902 if ( ( finehel->Prob() > 0.05 ) || nodrop )
903 {
904 MdcxTrklist.append( finehel );
905 // yzhang hist
906 for ( int iHit = 0; iHit < finehel->Nhits(); iHit++ )
907 {
908 if ( g_trklappend1 )
909 g_trklappend1->fill( finehel->Layer( iHit ) ); // yzhang hist
910 }
911 // zhangy hist
912 }
913 else
914 {
915 MdcxFittedHel* drophel = new MdcxFittedHel( drophits( finehel ) );
916 float ratdrop = float( drophel->Nhits() ) / float( finehel->Nhits() );
917 if ( 4 == m_debug )
918 cout << "APPEND track " << trkCount << ", drops "
919 << finehel->Nhits() - drophel->Nhits() << " helixnhit "
920 << finehel->Nhits() << " drophit " << drophel->Nhits()
921 << " hits, drop rate=" << ratdrop << ">?" << 0.5 << " prob "
922 << drophel->Prob() << " >?" << finehel->Prob() << " fail==0? "
923 << drophel->Fail() << endl;
924 if ( ( drophel->Prob() > finehel->Prob() ) && ( ratdrop > 0.50 ) &&
925 ( drophel->Fail() == 0 ) )
926 { // FIXME
927 if ( 4 == m_debug ) cout << "append drop helix " << endl;
928 MdcxTrklist.append( drophel );
929 if ( g_trklappend2 )
930 for ( int iHit = 0; iHit < drophel->Nhits(); iHit++ )
931 { g_trklappend2->fill( drophel->Layer( iHit ) ); }
932 delete finehel;
933 }
934 else
935 {
936 if ( 4 == m_debug ) cout << "append fine helix " << endl;
937 MdcxTrklist.append( finehel );
938 if ( g_trklappend3 )
939 for ( int iHit = 0; iHit < finehel->Nhits(); iHit++ )
940 { g_trklappend3->fill( finehel->Layer( iHit ) ); }
941 delete drophel;
942 }
943 }
944
945 // End hit dropping section
946 int nwl = MdcxTrklist.length() - 1;
947 if ( ( 4 == m_debug ) && ( nwl > -1 ) )
948 {
949 if ( 4 == m_debug )
950 {
951 cout << endl << "found track " << trkCount << std::endl;
952 MdcxTrklist[nwl]->print();
953 MdcxTrklist[nwl]->FitPrint();
954 }
955 }
956 }
957 else
958 { //*endif good fit so add to MdcxTrklist
959 for ( int i = 0; i < nSeg; i++ )
960 {
961 if ( segList[i]->GetUsedOnHel() == trkCount ) segList[i]->SetUsedOnHel( 0 );
962 }
963 delete finehel;
964 trkCount--;
965 }
966 } //*end if not duplicate
967 } //*end floop first_prob >= probmin
968 } // end iStV
969 } // end iStU
970 } // end iAx
971 if ( 4 == m_debug ) cout << "end of this combo" << endl;
972 } // end iSegCombo
973 if ( 4 == m_debug )
974 cout << " In MdcxFindTracks, found " << trkCount << " good tracks" << endl;
975} // end of process
AIDA::IHistogram1D * g_dPhiAV
NTuple::Item< double > m_addSegAddPhiLay
AIDA::IHistogram1D * g_trklappend2
NTuple::Item< double > m_addSegSeedPhi0
NTuple::Item< double > m_segCombDLenUV
AIDA::IHistogram1D * g_trklappend1
NTuple::Item< long > m_addSegSlayer
NTuple::Item< long > m_addSegSame
AIDA::IHistogram1D * g_trklproca
NTuple::Item< double > m_addSegSeedSl
NTuple::Tuple * m_xtupleAddSeg1
NTuple::Item< double > m_addSegLen
NTuple::Item< double > m_addSegSeedPhiLay
NTuple::Item< double > m_segCombSlV
NTuple::Item< double > m_addSegSeedD0
AIDA::IHistogram1D * g_trklappend3
AIDA::IHistogram1D * g_trklhelix
AIDA::IHistogram1D * g_trkllmk
NTuple::Item< double > m_segCombPhiA
NTuple::Item< double > m_segCombSlA
NTuple::Item< double > m_segCombDLenAU
AIDA::IHistogram1D * g_trklgood
NTuple::Item< double > m_segCombSlU
NTuple::Item< double > m_segCombPhiU
AIDA::IHistogram1D * g_dPhiAU
NTuple::Item< double > m_segCombSameUV
NTuple::Item< double > m_addSegAddPhi
NTuple::Item< double > m_addSegSeedPhi
NTuple::Item< double > m_segCombPhiV
NTuple::Item< double > m_addSegAddPhi0
NTuple::Tuple * m_xtupleSegComb
NTuple::Item< double > m_addSegAddSl
NTuple::Item< double > m_segCombOmega
NTuple::Item< double > m_addSegAddD0
NTuple::Item< double > m_segCombSameAU
AIDA::IHistogram1D * g_trklfirstProb
NTuple::Tuple * m_xtupleAddSeg2
AIDA::IHistogram1D * g_omegag
NTuple::Item< long > m_addSegEvtNo
NTuple::Item< double > m_addSegPoca
float Mdcxprobab(int &ndof, float &chisq)
Definition Mdcxprobab.cxx:4
#define prt(n)
MdcxFittedHel drophits(MdcxFittedHel *f)
double findz_sl(int i1, int i2, const HepAList< MdcxSeg > &slclus)
double findz_cyl(int i1, int i2, const HepAList< MdcxSeg > &slclus)
bool testFromSameTrack(MdcxSeg *seg1, MdcxSeg *seg2)
double dlen(int slay1, double p1, int slay2, double p2, double om)
MdcxHel TakeToOrigin(MdcxHel &)
int Layer(int hitno=0) const
int Fail(float Probmin=0.0) const
MdcxFittedHel & Grow(const MdcxFittedHel &, const HepAList< MdcxHit > &)
double X0() const
Definition MdcxHel.cxx:85
double Lmax() const
Definition MdcxHel.cxx:137
double Yc() const
Definition MdcxHel.cxx:76
double Xc() const
Definition MdcxHel.cxx:67
void print() const
Definition MdcxHel.cxx:403
double Y0() const
Definition MdcxHel.cxx:87
void flip()
Definition MdcxHel.cxx:421
static const int nSegCombo
relative to MdcxFindTracks
static const int findTrkGroup[nSegCombo][3]
– relative to MdcxFindTracks

Referenced by MdcxFindTracks().

◆ process() [2/3]

void MdcxFindTracks::process ( const HepAList< MdcxSeg > & f)

◆ process() [3/3]

void MdcxFindTracks::process ( const HepAList< MdcxSeg > & f)

◆ resout() [1/3]

void MdcxFindTracks::resout ( MdcxFittedHel * f)
protected

Definition at line 1123 of file MdcxFindTracks.cxx.

1123 {
1124 MdcxHel mdcxHel = (MdcxHel)*finehel;
1125 HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
1126 int kkk = 0;
1127 while ( hitlist[kkk] )
1128 {
1129 hitlist[kkk]->SetConstErr( 0 );
1130 kkk++;
1131 }
1132 MdcxFittedHel newhel( hitlist, mdcxHel );
1133 if ( 4 == m_debug ) newhel.FitPrint( mhel, cout );
1134 kkk = 0;
1135 while ( hitlist[kkk] )
1136 {
1137 float doca = newhel.Doca( *hitlist[kkk] );
1138 float zh = newhel.Doca_Zh();
1139 float tof = newhel.Doca_Tof();
1140 float dd = hitlist[kkk]->d( newhel );
1141 float tt = hitlist[kkk]->tcor( zh, tof );
1142 if ( 4 == m_debug )
1143 {
1144 cout << "MdcxFindTracks::resout (" << hitlist[kkk]->Layer() << ","
1145 << hitlist[kkk]->WireNo() << ") dt " << tt << " resid " << fabs( dd ) - fabs( doca )
1146 << " doca " << fabs( doca ) << endl;
1147 }
1148 kkk++;
1149 }
1150 kkk = 0;
1151 while ( hitlist[kkk] )
1152 {
1153 hitlist[kkk]->SetConstErr( 1 );
1154 kkk++;
1155 }
1156} // endof drophits

◆ resout() [2/3]

void MdcxFindTracks::resout ( MdcxFittedHel * f)
protected

◆ resout() [3/3]

void MdcxFindTracks::resout ( MdcxFittedHel * f)
protected

◆ TakeToOrigin() [1/3]

MdcxHel MdcxFindTracks::TakeToOrigin ( MdcxHel & ihel)
protected

Definition at line 1158 of file MdcxFindTracks.cxx.

1158 {
1159 double lng;
1160 double omegag = ihel.Omega();
1161 double phi0g = ihel.Phi0();
1162 double phi0s = phi0g;
1163 double xp = ihel.Xref() - ihel.SinPhi0() * ihel.D0();
1164 double yp = ihel.Yref() + ihel.CosPhi0() * ihel.D0();
1165 double d0g = yp * cos( phi0g ) - xp * sin( phi0g );
1166 if ( omegag != 0.0 )
1167 {
1168 phi0g = phi0g - Constants::pi;
1169 double xc = sin( phi0g ) / omegag;
1170 double yc = -cos( phi0g ) / omegag;
1171 double t1 = -xc - xp;
1172 double t2 = yc + yp;
1173 double phinew = atan2( t1, t2 );
1174 if ( omegag > 0.0 ) phinew += Constants::pi;
1175 if ( phinew > Constants::pi ) phinew -= Constants::twoPi;
1176 double xh = xc - sin( phinew ) / omegag;
1177 double yh = yc + cos( phinew ) / omegag;
1178 double x0 = xh + xp;
1179 double y0 = yh + yp;
1180 d0g = sqrt( x0 * x0 + y0 * y0 );
1181 phi0g = phinew + Constants::pi; ////???
1182 double sd0 = x0 * sin( phi0g ) - y0 * cos( phi0g );
1183 if ( sd0 > 0.0 ) d0g = -d0g;
1184 lng = dlen( -2, phi0g, -1, phi0s, omegag ); // yzhang TEMP 2011-10-17
1185 // lng = dlen(phi0g, phi0s, omegag);
1186 }
1187 else { lng = ihel.Xref() * ihel.CosPhi0() + ihel.Yref() * ihel.SinPhi0(); }
1188 double tanlg = ihel.Tanl();
1189 double z0g = ihel.Z0() - lng * tanlg;
1190 // cout << " z0g, tanlg, lng " << z0g << " " << tanlg << " " << lng << endl;
1191 double t0g = ihel.T0();
1192 int codeg = ihel.Code();
1193 MdcxHel ohel( d0g, phi0g, omegag, z0g, tanlg, t0g, codeg );
1194 return ohel;
1195}

Referenced by process().

◆ TakeToOrigin() [2/3]

MdcxHel MdcxFindTracks::TakeToOrigin ( MdcxHel & )
protected

◆ TakeToOrigin() [3/3]

MdcxHel MdcxFindTracks::TakeToOrigin ( MdcxHel & )
protected

◆ testFromSameTrack() [1/3]

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

Definition at line 1197 of file MdcxFindTracks.cxx.

1197 {
1198 /*
1199 std::cout<< " tkid seg1 "<<std::endl;
1200 for (int i =0; i<seg1->Nhits(); i++){
1201 int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
1202 std::cout<< tkid <<std::endl;
1203 }
1204 std::cout<< " tkid seg2 "<<std::endl;
1205 for (int i =0; i<seg2->Nhits(); i++){
1206 int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
1207 std::cout<< tkid <<std::endl;
1208 }
1209 */
1210
1211 int trackId = -9999;
1212 for ( int i = 0; i < seg1->Nhits(); i++ )
1213 {
1214 int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
1215 if ( tkid >= 1000 ) tkid -= 1000;
1216 if ( ( trackId == -9999 ) )
1217 {
1218 if ( tkid >= 0 ) trackId = tkid;
1219 }
1220 else
1221 {
1222 if ( tkid != trackId ) return false;
1223 }
1224 }
1225 for ( int i = 0; i < seg2->Nhits(); i++ )
1226 {
1227 int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
1228 if ( tkid >= 1000 ) tkid -= 1000;
1229 if ( ( trackId == -9999 ) ) { return false; }
1230 else
1231 {
1232 if ( tkid != trackId ) return false;
1233 }
1234 }
1235 return true;
1236}

Referenced by process().

◆ testFromSameTrack() [2/3]

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

◆ testFromSameTrack() [3/3]

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

Member Data Documentation

◆ curvmax

float MdcxFindTracks::curvmax
protected

◆ m_debug

int MdcxFindTracks::m_debug
protected

◆ m_doSag

bool MdcxFindTracks::m_doSag
protected

◆ MdcxTrklist

HepAList< MdcxFittedHel > MdcxFindTracks::MdcxTrklist
protected

◆ mhel

MdcxHel MdcxFindTracks::mhel
protected

◆ nbad

int MdcxFindTracks::nbad
protected

◆ nfail

int MdcxFindTracks::nfail
protected

◆ ngood

int MdcxFindTracks::ngood
protected

◆ nhitsgd

int MdcxFindTracks::nhitsgd
protected

◆ nhitsmc

int MdcxFindTracks::nhitsmc
protected

◆ npure

int MdcxFindTracks::npure
protected

◆ probmin

float MdcxFindTracks::probmin
protected

The documentation for this class was generated from the following files: