BOSS 6.6.4.p01
BESIII Offline Software System
Loading...
Searching...
No Matches
Public Member Functions | Protected Member Functions | Protected Attributes | List of all members
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
 

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() [1/2]

MdcxFindTracks::MdcxFindTracks ( )

Definition at line 110 of file MdcxFindTracks.cxx.

110 {
111 probmin = MdcxParameters::minTrkProb;
112 curvmax = 1000000.0;
113}
MdcxParameters::minTrkProb
static double minTrkProb
Definition: MdcxParameters.h:15

◆ MdcxFindTracks() [2/2]

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

Definition at line 115 of file MdcxFindTracks.cxx.

115 {
116 probmin = MdcxParameters::minTrkProb;
117 curvmax = 1000000.0;
118 m_debug = debug;
119 process(segList);
120}
MdcxFindTracks::process
void process(const HepAList< MdcxSeg > &f)
Definition: MdcxFindTracks.cxx:124

◆ ~MdcxFindTracks()

MdcxFindTracks::~MdcxFindTracks ( )
virtual

Definition at line 122 of file MdcxFindTracks.cxx.

122{ KillList(); }

Member Function Documentation

◆ beginmore()

void MdcxFindTracks::beginmore ( )
protected

Definition at line 839 of file MdcxFindTracks.cxx.

839 {
840 ngood=0; nbad=0; nfail=0; npure=0; // zero scalars
841 nhitsmc=0; nhitsgd=0; // zero scalars
842 printpar(cout); // print parameters to screen
843}//endof beginmore
MdcxFindTracks::printpar
void printpar(std::ostream &o)
Definition: MdcxFindTracks.cxx:833

◆ dlen()

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

Definition at line 894 of file MdcxFindTracks.cxx.

894 {
895 //double dp = p2 - p1;
896 double dp;//yzhang change 2011-10-17
897 if((slay1==2 || slay1==3 || slay1 == 4 || slay1 == 9 || slay1 ==10) ){ //Ax
898 dp = p2 - p1;
899 }else if((slay2 >= slay1)) {
900 dp = p2 - p1;
901 }else {
902 dp = p1 - p2;
903 }
904
905 if (om < 0.0) {
906 while (dp < -Constants::pi) dp += Constants::pi;
907 while (dp > 0.0) dp -= 1*Constants::pi;
908 } else {
909 while (dp < 0.0) dp += 1*Constants::pi;
910 while (dp > Constants::pi) dp -= 1*Constants::pi;
911 }
912 //double dl = dp * r;
913 double dl = dp/om;
914 //if(4 == m_debug) {
915 // double r = 1./om;
916 // cout prt(p1) prt(p2) prt(dp) prt(om) prt(r) prt(dl) << endl;
917 //}
918 return dl;
919}
Constants::pi
static const double pi
Definition: Constants.h:38

Referenced by process(), and TakeToOrigin().

◆ drophits()

MdcxFittedHel MdcxFindTracks::drophits ( MdcxFittedHel *  f)
protected

Definition at line 921 of file MdcxFindTracks.cxx.

921 {
922 MdcxHel mdcxHel = (MdcxHel) *finehel;
923 const HepAList<MdcxHit>& hitlist = finehel->XHitList();
924 HepAList<MdcxHit> nwhitlist = hitlist;
925 int nremove = 0;
926 int kkk = 0;
927 while(nwhitlist[kkk]) {
928 double pull = nwhitlist[kkk]->pull(mdcxHel);
929 int layer = nwhitlist[kkk]->Layer();
930
931 if(m_xtupleDropHits){
932 float t_doca = mdcxHel.Doca(*(nwhitlist[kkk]));
933 float t_e = nwhitlist[kkk]->e(t_doca);
934 m_segDropHitsEvtNo = g_eventNo;
935 m_segDropHitsLayer = layer;
936 m_segDropHitsWire = nwhitlist[kkk]->WireNo();
937 m_segDropHitsPull = pull;
938 m_segDropHitsDoca = t_doca;
939 m_segDropHitsSigma = t_e;
940 m_segDropHitsDrift = nwhitlist[kkk]->d(mdcxHel);
941 m_segDropHitsMcTkId = nwhitlist[kkk]->getDigi()->getTrackIndex();
942 m_xtupleDropHits->write();
943 }
944
945 if (fabs(pull) > MdcxParameters::dropHitsSigma[layer]){
946 nwhitlist.remove(kkk);
947 nremove++;
948 } else {
949 kkk++;
950 }
951 }
952 if(m_debug==4) cout<< " MdcxFindTracks::drophits drop "<<nremove<<" hits "
953 << " nhits of drop helix = "<<nwhitlist.length()<<endl;
954
955 MdcxFittedHel newhel(nwhitlist,mdcxHel);
956 return newhel;
957}//endof drophits
m_segDropHitsWire
NTuple::Item< long > m_segDropHitsWire
Definition: MdcxHistItem.h:194
m_segDropHitsDoca
NTuple::Item< double > m_segDropHitsDoca
Definition: MdcxHistItem.h:196
m_segDropHitsSigma
NTuple::Item< double > m_segDropHitsSigma
Definition: MdcxHistItem.h:197
m_segDropHitsLayer
NTuple::Item< long > m_segDropHitsLayer
Definition: MdcxHistItem.h:193
m_segDropHitsMcTkId
NTuple::Item< double > m_segDropHitsMcTkId
Definition: MdcxHistItem.h:199
m_segDropHitsEvtNo
NTuple::Item< long > m_segDropHitsEvtNo
Definition: MdcxHistItem.h:192
m_segDropHitsPull
NTuple::Item< double > m_segDropHitsPull
Definition: MdcxHistItem.h:195
m_segDropHitsDrift
NTuple::Item< double > m_segDropHitsDrift
Definition: MdcxHistItem.h:198
m_xtupleDropHits
NTuple::Tuple * m_xtupleDropHits
Definition: MdcxHistItem.h:191
MdcxHel::Doca
double Doca(double WX, double WY, double WZ, double X, double Y, double Z=0.0)
Definition: MdcxHel.cxx:239
MdcxParameters::dropHitsSigma
static float dropHitsSigma[43]
Definition: MdcxParameters.h:43
g_eventNo
NTuple::Item< long > g_eventNo
Definition: ntupleItem.h:22

Referenced by process().

◆ endmore()

void MdcxFindTracks::endmore ( )
protected

Definition at line 845 of file MdcxFindTracks.cxx.

845 {
846 print(cout);
847}//endof endmore
MdcxFindTracks::print
void print(std::ostream &o)
Definition: MdcxFindTracks.cxx:821

◆ findz_cyl()

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

Definition at line 869 of file MdcxFindTracks.cxx.

869 {
870 double zint = -1000.0;
871 double r = 1.0/fabs(segList[iAx]->Omega());
872 double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0();
873 double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
874 double sx = segList[iStU]->Xline_slope();
875 double sy = segList[iStU]->Yline_slope();
876 double cx = xl_0 - segList[iAx]->Xc();
877 double cy = yl_0 - segList[iAx]->Yc();
878 double a = sx*sx + sy*sy;
879 double b = 2.0*(sx*cx + sy*cy);
880 //double c = (cx*cx + cy*cy - r*r);
881 double bsq = b * b;
882 double fac = 4.0 * a * (cx*cx + cy*cy - r*r);
883 double ta = 2.0 * a;
884 if (fac < bsq) {
885 double sfac = sqrt(bsq - fac);
886 double z1 = -(b - sfac)/ta;
887 double z2 = -(b + sfac)/ta;
888 zint = (fabs(z1) < fabs(z2)) ? z1 : z2;
889 }
890 // cout << " findz_cyl segs " << iAx << " " << iStU << " zint " << zint << endl;
891 return zint;
892}

Referenced by process().

◆ findz_sl()

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

Definition at line 849 of file MdcxFindTracks.cxx.

849 {
850 double Ap = -sin(segList[iAx]->Phi0_sl_approx());
851 double Bp = cos(segList[iAx]->Phi0_sl_approx());
852 double xp = segList[iAx]->Xref() + Ap*segList[iAx]->D0_sl_approx();
853 double yp = segList[iAx]->Yref() + Bp*segList[iAx]->D0_sl_approx();
854 double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0();
855 double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
856
857 const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
858 double Cl = hitsStU[0]->wz();
859 double ndotl = Bp*hitsStU[0]->wy() + Ap*hitsStU[0]->wx();
860 /*
861 std::cout<<"ndotl "<<ndotl<<" "<<Bp<<" "<<hitsStU[0]->wy()<<" "<<
862 Ap<<" "<<hitsStU[0]->wx()<<std::endl;
863 */
864 double zint = (ndotl != 0.0) ? ((Bp*(yp-yl)+Ap*(xp-xl))*Cl/ndotl) : -1000.;
865 // cout << " findz_sl segs " << iAx << " " << iStU << " zint " << zint << endl;
866 return zint;
867}
sin
double sin(const BesAngle a)
Definition: BesAngle.h:210
cos
double cos(const BesAngle a)
Definition: BesAngle.h:213

Referenced by process().

◆ GetMdcxTrklist()

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

Definition at line 39 of file MdcxFindTracks.h.

39{return MdcxTrklist;}
MdcxFindTracks::MdcxTrklist
HepAList< MdcxFittedHel > MdcxTrklist
Definition: MdcxFindTracks.h:47

Referenced by MdcxTrackFinder::execute().

◆ KillList()

void MdcxFindTracks::KillList ( )
inlineprotected

Definition at line 61 of file MdcxFindTracks.h.

61{HepAListDeleteAll(MdcxTrklist);}

Referenced by ~MdcxFindTracks().

◆ plot()

void MdcxFindTracks::plot ( ) const

Definition at line 830 of file MdcxFindTracks.cxx.

830 {
831} // endof plot

◆ print()

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

Definition at line 821 of file MdcxFindTracks.cxx.

821 {
822 o<<"ngood:"<<ngood<<endl;
823 o<<"npure:"<<npure<<endl;
824 o<<"nbad:"<<nbad<<endl;
825 o<<"nfail:"<<nfail<<endl;
826 o<<"nhitsmc:"<<nhitsmc<<endl;
827 o<<"nhitsgd:"<<nhitsgd<<endl;
828} // endof print

Referenced by endmore().

◆ printpar()

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

Definition at line 833 of file MdcxFindTracks.cxx.

833 {
834 o << "|MdcxFindTracks| parameters" << endl;
835 o << "probmin:" << probmin << endl;
836 o << "curvmax:" << curvmax << endl;
837}//endof printpar

Referenced by beginmore().

◆ process()

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

Referenced by MdcxFindTracks().

◆ resout()

void MdcxFindTracks::resout ( MdcxFittedHel *  f)
protected

Definition at line 959 of file MdcxFindTracks.cxx.

959 {
960 MdcxHel mdcxHel = (MdcxHel) *finehel;
961 HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
962 int kkk = 0; while(hitlist[kkk]){hitlist[kkk]->SetConstErr(0); kkk++;}
963 MdcxFittedHel newhel(hitlist, mdcxHel);
964 if(4 == m_debug) newhel.FitPrint(mhel,cout);
965 kkk = 0;
966 while(hitlist[kkk]) {
967 float doca = newhel.Doca(*hitlist[kkk]);
968 float zh = newhel.Doca_Zh();
969 float tof = newhel.Doca_Tof();
970 float dd = hitlist[kkk]->d(newhel);
971 float tt = hitlist[kkk]->tcor(zh, tof);
972 if (4 == m_debug){
973 cout << "MdcxFindTracks::resout ("
974 << hitlist[kkk]->Layer() << ","
975 << hitlist[kkk]->WireNo() << ") dt "
976 << tt << " resid "
977 << fabs(dd)-fabs(doca) << " doca "
978 << fabs(doca) << endl;
979 }
980 kkk++;
981 }
982 kkk=0;
983 while(hitlist[kkk]){hitlist[kkk]->SetConstErr(1); kkk++;}
984}//endof drophits

◆ TakeToOrigin()

MdcxHel MdcxFindTracks::TakeToOrigin ( MdcxHel &  ihel)
protected

Definition at line 986 of file MdcxFindTracks.cxx.

986 {
987 double lng;
988 double omegag = ihel.Omega();
989 double phi0g = ihel.Phi0();
990 double phi0s = phi0g;
991 double xp = ihel.Xref() - ihel.SinPhi0()*ihel.D0();
992 double yp = ihel.Yref() + ihel.CosPhi0()*ihel.D0();
993 double d0g = yp*cos(phi0g) - xp*sin(phi0g);
994 if (omegag != 0.0) {
995 phi0g = phi0g - Constants::pi;
996 double xc = sin(phi0g)/omegag;
997 double yc = -cos(phi0g)/omegag;
998 double t1 = -xc - xp;
999 double t2 = yc + yp;
1000 double phinew = atan2(t1, t2);
1001 if (omegag > 0.0) phinew += Constants::pi;
1002 if (phinew > Constants::pi) phinew -= Constants::twoPi;
1003 double xh = xc - sin(phinew)/omegag;
1004 double yh = yc + cos(phinew)/omegag;
1005 double x0 = xh + xp;
1006 double y0 = yh + yp;
1007 d0g = sqrt(x0*x0 + y0*y0);
1008 phi0g = phinew + Constants::pi; ////???
1009 double sd0 = x0*sin(phi0g) - y0*cos(phi0g);
1010 if (sd0 > 0.0) d0g = -d0g;
1011 lng = dlen(-2, phi0g, -1, phi0s, omegag);//yzhang TEMP 2011-10-17
1012 //lng = dlen(phi0g, phi0s, omegag);
1013 } else {
1014 lng = ihel.Xref()*ihel.CosPhi0() + ihel.Yref()*ihel.SinPhi0();
1015 }
1016 double tanlg = ihel.Tanl();
1017 double z0g = ihel.Z0() - lng*tanlg;
1018 //cout << " z0g, tanlg, lng " << z0g << " " << tanlg << " " << lng << endl;
1019 double t0g = ihel.T0();
1020 int codeg = ihel.Code();
1021 MdcxHel ohel(d0g, phi0g, omegag, z0g, tanlg, t0g, codeg);
1022 return ohel;
1023}
MdcxHel::CosPhi0
double CosPhi0() const
Definition: MdcxHel.h:63
MdcxHel::T0
double T0() const
Definition: MdcxHel.h:62
MdcxHel::Code
int Code() const
Definition: MdcxHel.h:74
MdcxHel::Yref
double Yref() const
Definition: MdcxHel.h:61
MdcxHel::Z0
double Z0() const
Definition: MdcxHel.h:56
MdcxHel::Tanl
double Tanl() const
Definition: MdcxHel.h:57
MdcxHel::SinPhi0
double SinPhi0() const
Definition: MdcxHel.h:64
MdcxHel::Xref
double Xref() const
Definition: MdcxHel.h:59

Referenced by process().

◆ testFromSameTrack()

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

Definition at line 1025 of file MdcxFindTracks.cxx.

1025 {
1026 /*
1027 std::cout<< " tkid seg1 "<<std::endl;
1028 for (int i =0; i<seg1->Nhits(); i++){
1029 int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
1030 std::cout<< tkid <<std::endl;
1031 }
1032 std::cout<< " tkid seg2 "<<std::endl;
1033 for (int i =0; i<seg2->Nhits(); i++){
1034 int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
1035 std::cout<< tkid <<std::endl;
1036 }
1037 */
1038
1039 int trackId = -9999;
1040 for (int i =0; i<seg1->Nhits(); i++){
1041 int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
1042 if(tkid>=1000) tkid -=1000;
1043 if((trackId == -9999)){
1044 if(tkid>=0) trackId = tkid;
1045 }else{
1046 if(tkid!=trackId) return false;
1047 }
1048 }
1049 for (int i =0; i<seg2->Nhits(); i++){
1050 int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
1051 if(tkid>=1000) tkid -=1000;
1052 if((trackId == -9999)){
1053 return false;
1054 }else{
1055 if(tkid!=trackId) return false;
1056 }
1057 }
1058 return true;
1059}

Referenced by process().

Member Data Documentation

◆ curvmax

float MdcxFindTracks::curvmax
protected

Definition at line 50 of file MdcxFindTracks.h.

Referenced by MdcxFindTracks(), and printpar().

◆ m_debug

int MdcxFindTracks::m_debug
protected

Definition at line 57 of file MdcxFindTracks.h.

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

◆ m_doSag

bool MdcxFindTracks::m_doSag
protected

Definition at line 58 of file MdcxFindTracks.h.

◆ MdcxTrklist

HepAList<MdcxFittedHel> MdcxFindTracks::MdcxTrklist
protected

Definition at line 47 of file MdcxFindTracks.h.

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

◆ mhel

MdcxHel MdcxFindTracks::mhel
protected

Definition at line 48 of file MdcxFindTracks.h.

Referenced by process(), and resout().

◆ nbad

int MdcxFindTracks::nbad
protected

Definition at line 53 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ nfail

int MdcxFindTracks::nfail
protected

Definition at line 54 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ ngood

int MdcxFindTracks::ngood
protected

Definition at line 51 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ nhitsgd

int MdcxFindTracks::nhitsgd
protected

Definition at line 56 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ nhitsmc

int MdcxFindTracks::nhitsmc
protected

Definition at line 55 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ npure

int MdcxFindTracks::npure
protected

Definition at line 52 of file MdcxFindTracks.h.

Referenced by beginmore(), and print().

◆ probmin

float MdcxFindTracks::probmin
protected

Definition at line 49 of file MdcxFindTracks.h.

Referenced by MdcxFindTracks(), printpar(), and process().

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