00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include <math.h>
00023 #include "MdcxReco/MdcxFindTracks.h"
00024 #include "MdcxReco/MdcxParameters.h"
00025 #include "MdcxReco/Mdcxprobab.h"
00026 #include "MdcxReco/MdcxHit.h"
00027 #include "MdcGeom/Constants.h"
00028
00029 #include "AIDA/IHistogram1D.h"
00030 #include "AIDA/IHistogram2D.h"
00031 #include "GaudiKernel/NTuple.h"
00032
00033 #define GET_NAME(n) #n
00034 #define prt(n) <<setw(8)<<GET_NAME(n)<<" "<<setw(8)<<n<<endl
00035
00036 extern int g_eventNo;
00037 extern AIDA::IHistogram1D* g_omegag;
00038 extern AIDA::IHistogram2D* g_dPhiAU;
00039 extern AIDA::IHistogram2D* g_dPhiAV;
00040 extern AIDA::IHistogram1D* g_trkllmk;
00041 extern AIDA::IHistogram1D* g_trklcircle;
00042 extern AIDA::IHistogram1D* g_trklgood;
00043 extern AIDA::IHistogram1D* g_trklhelix;
00044 extern AIDA::IHistogram1D* g_trkldrop1;
00045 extern AIDA::IHistogram1D* g_trkldrop2;
00046 extern AIDA::IHistogram1D* g_trklappend1;
00047 extern AIDA::IHistogram1D* g_trklappend2;
00048 extern AIDA::IHistogram1D* g_trklappend3;
00049 extern AIDA::IHistogram1D* g_trklfirstProb;
00050 extern AIDA::IHistogram1D* g_trkltemp;
00051 extern AIDA::IHistogram1D* g_trklproca;
00052 extern AIDA::IHistogram1D* g_trkld;
00053 extern AIDA::IHistogram1D* g_trkle;
00054 extern AIDA::IHistogram1D* g_trkldoca;
00055 extern AIDA::IHistogram1D* g_trkllayer;
00056 extern AIDA::IHistogram2D* g_trkldl;
00057 extern AIDA::IHistogram2D* g_trklel;
00058 extern AIDA::IHistogram2D* g_dropHitsSigma;
00059
00060 extern NTuple::Tuple* m_xtupleAddSeg1;
00061 extern NTuple::Item<long> m_addSegSame;
00062 extern NTuple::Item<double> m_addSegSeedSl;
00063 extern NTuple::Item<double> m_addSegSeedPhi;
00064 extern NTuple::Item<double> m_addSegSeedPhiLay;
00065 extern NTuple::Item<double> m_addSegSeedD0;
00066 extern NTuple::Item<double> m_addSegSeedLen;
00067 extern NTuple::Item<double> m_addSegSeedPhi0;
00068 extern NTuple::Item<double> m_addSegAddSl;
00069 extern NTuple::Item<double> m_addSegAddPhi;
00070 extern NTuple::Item<double> m_addSegAddPhiLay;
00071 extern NTuple::Item<double> m_addSegAddD0;
00072 extern NTuple::Item<double> m_addSegAddLen;
00073 extern NTuple::Item<double> m_addSegAddPhi0;
00074
00075 extern NTuple::Tuple* m_xtupleAddSeg2;
00076 extern NTuple::Item<long> m_addSegEvtNo;
00077 extern NTuple::Item<double> m_addSegPoca;
00078 extern NTuple::Item<long> m_addSegSlayer;
00079 extern NTuple::Item<double> m_addSegLen;
00080
00081 extern NTuple::Tuple* m_xtupleSegComb;
00082 extern NTuple::Item<double> m_segCombOmega;
00083 extern NTuple::Item<double> m_segCombSameAU;
00084 extern NTuple::Item<double> m_segCombSameUV;
00085 extern NTuple::Item<double> m_segCombDLenAU;
00086 extern NTuple::Item<double> m_segCombDLenUV;
00087 extern NTuple::Item<double> m_segCombSlA;
00088 extern NTuple::Item<double> m_segCombSlU;
00089 extern NTuple::Item<double> m_segCombSlV;
00090 extern NTuple::Item<double> m_segCombPhiA;
00091 extern NTuple::Item<double> m_segCombPhiU;
00092 extern NTuple::Item<double> m_segCombPhiV;
00093 extern NTuple::Item<double> m_segCombPoca;
00094
00095 extern NTuple::Tuple* m_xtupleDropHits;
00096 extern NTuple::Item<long> m_segDropHitsEvtNo;
00097 extern NTuple::Item<long> m_segDropHitsLayer;
00098 extern NTuple::Item<long> m_segDropHitsWire;
00099 extern NTuple::Item<double> m_segDropHitsPull;
00100 extern NTuple::Item<double> m_segDropHitsDoca;
00101 extern NTuple::Item<double> m_segDropHitsSigma;
00102 extern NTuple::Item<double> m_segDropHitsDrift;
00103 extern NTuple::Item<double> m_segDropHitsMcTkId;
00104
00105
00106 using std::cout;
00107 using std::endl;
00108 using std::ostream;
00109
00110 MdcxFindTracks::MdcxFindTracks(){
00111 probmin = MdcxParameters::minTrkProb;
00112 curvmax = 1000000.0;
00113 }
00114
00115 MdcxFindTracks::MdcxFindTracks(const HepAList<MdcxSeg> &segList, int debug){
00116 probmin = MdcxParameters::minTrkProb;
00117 curvmax = 1000000.0;
00118 m_debug = debug;
00119 process(segList);
00120 }
00121
00122 MdcxFindTracks::~MdcxFindTracks() { KillList(); }
00123
00124 void MdcxFindTracks::process(const HepAList<MdcxSeg> &segList) {
00125 static int nfound = 0;
00126 static float sumprob = 0.0;
00127
00128 MdcxHel thel(0., 0., 0., 0., 0.);
00129 mhel = thel;
00130 int nSeg = segList.length();
00131 for (int i = 0; i < nSeg ; i++) { segList[i]->SetUsedOnHel(0); }
00132
00133 double xp, yp, xl1, yl1, rl1, xl2, yl2, rl2;
00134 double d0g, phi0g, omegag, z0g, tanlg;
00135 double d0g_sl, phi0g_sl, omegag_sl, z0g_sl, tanlg_sl;
00136 double zintAU, zintAV, zintAU_sl, zintAV_sl = 9999.;
00137 double rl1_sl, rl2_sl;
00138 double xc, yc;
00139 double sinPhi0g_sl, cosPhi0g_sl, xp_sl, yp_sl;
00140 double phiAx, phiStU, phiStV, xl1_sl, yl1_sl, xl2_sl, yl2_sl, xl1_0, yl1_0, xl2_0, yl2_0;
00141 double sx1, sy1, sx2, sy2;
00142 double x0g, y0g;
00143 double fltAx, ztest, fltL_sl, ztest_sl;
00144 int floop, trkCount= 0;
00145
00146
00147 for (int iSegCombo = 0; iSegCombo < MdcxParameters::nSegCombo; iSegCombo++) {
00148
00149 int axlay = MdcxParameters::findTrkGroup[iSegCombo][0];
00150 int stUlay = MdcxParameters::findTrkGroup[iSegCombo][1];
00151 int stVlay = MdcxParameters::findTrkGroup[iSegCombo][2];
00152
00153 if(4 == m_debug) std::cout <<" Test combo: ("<< axlay<<","<<stUlay<<","<<stVlay <<")"<< std::endl;
00154
00155 float max_dPhiAU = MdcxParameters::maxDp12[iSegCombo];
00156 float max_dPhiAV = MdcxParameters::maxDp13[iSegCombo];
00157
00158
00159
00160
00161 for (int iAx = 0; iAx < nSeg; iAx++) {
00162 bool b_useGood_stU = true;
00163 bool b_useGood_stV = true;
00164 if ((segList[iAx]->GetUsedOnHel() != 0) || (segList[iAx]->Fail() != 0)
00165 || (segList[iAx]->SuperLayer(0) != axlay) || (segList[iAx]->Origin() != -1) ) continue;
00166
00167 if (4 == m_debug){
00168 std::cout<< "---------1.Ax seg------ No."<<iAx<< std::endl;
00169 segList[iAx]->printSegAll();
00170 }
00171
00172
00173 omegag = segList[iAx]->Omega();
00174 if (fabs(omegag) > MdcxParameters::maxTrkOmega) {
00175 if (4 == m_debug) std::cout<< "SKIP by maxTrkOmega "<<fabs(omegag) << std::endl;
00176 continue;
00177 }
00178 if(4 == m_debug) std::cout << " iAx omegag = " << omegag << std::endl;
00179 if(g_omegag) g_omegag->fill(omegag);
00180
00181 const HepAList<MdcxHit>& hitsAx = segList[iAx]->XHitList();
00182
00183
00184 d0g_sl = segList[iAx]->D0_sl_approx();
00185 phi0g_sl = segList[iAx]->Phi0_sl_approx();
00186 omegag_sl = 0.0;
00187
00188
00189 sinPhi0g_sl = -sin(phi0g_sl);
00190 cosPhi0g_sl = cos(phi0g_sl);
00191 xp_sl = segList[iAx]->Xref() + sinPhi0g_sl*d0g_sl;
00192 yp_sl = segList[iAx]->Yref() + cosPhi0g_sl*d0g_sl;
00193 d0g_sl = yp_sl*cosPhi0g_sl + xp_sl*sinPhi0g_sl;
00194 fltL_sl = segList[iAx]->Xref()*cosPhi0g_sl - segList[iAx]->Yref()*sinPhi0g_sl;
00195 if(4 == m_debug) {
00196 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;
00197 }
00198
00199 if (fabs(omegag) > Constants::epsilon) {
00200 MdcxHel ohel = TakeToOrigin(*segList[iAx]);
00201 omegag = ohel.Omega();
00202 phi0g = ohel.Phi0();
00203 d0g = ohel.D0();
00204 xc = ohel.Xc();
00205 yc = ohel.Yc();
00206 x0g = ohel.X0();
00207 y0g = ohel.Y0();
00208 phiAx = atan2(y0g-yc, x0g-xc);
00209 fltAx = dlen(-2, phi0g,-1, segList[iAx]->Phi0(), omegag);
00210
00211 if(4 == m_debug) {
00212 std::cout <<"--Ax :"<< " ohel ";
00213 ohel.print();
00214 }
00215 }
00216
00217
00218
00219
00220
00221 for (int iStU = -1; true; ) {
00222 if(!b_useGood_stU && (iStU >= (nSeg-1))) break;
00223 if(b_useGood_stU && (iStU >= (nSeg-1))){
00224 iStU = -1;
00225 b_useGood_stU = false;
00226 }
00227 iStU++;
00228 if ((segList[iAx]->GetUsedOnHel() != 0)
00229 || (segList[iStU]->GetUsedOnHel() != 0)
00230 || (segList[iStU]->Fail() != 0)
00231 || (segList[iStU]->SuperLayer(0) != stUlay)
00232 || (segList[iStU]->Origin() != -1)
00233
00234 || (b_useGood_stU && (segList[iStU]->Nhits()< 4))
00235 ){ continue;}
00236
00237 if (4 == m_debug){
00238 std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl;
00239 std::cout<< "---------2.St U seg ------No."<<iStU<< std::endl;
00240 segList[iStU]->printSeg();
00241 std::cout<<" omega of seg iStU = "<<segList[iStU]->Omega()<<std::endl;
00242 }
00243
00244
00245 const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
00246 double dPhiAU = fabs(hitsAx[0]->phiMid() - hitsStU[0]->phiMid());
00247 if (dPhiAU > Constants::pi) dPhiAU = Constants::twoPi - dPhiAU;
00248 if(4 == m_debug) {
00249 if (dPhiAU > max_dPhiAU){
00250 cout << "SKIP by dPhiAU " <<dPhiAU <<" > "<< max_dPhiAU << endl;
00251 }else{
00252 cout << "KEEP by dPhiAU " <<dPhiAU<< " <= " << max_dPhiAU << endl;
00253 }
00254 std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
00255 <<" StU "<<hitsStU[0]->phiMid() << std::endl;
00256 std::cout<< std::endl;
00257 }
00258 if (g_dPhiAU) g_dPhiAU->fill(dPhiAU,segList[iStU]->SuperLayer());
00259
00260 if (dPhiAU > max_dPhiAU){ continue; }
00261
00262
00263 xl1_0 = segList[iStU]->Xline_bbrrf();
00264 yl1_0 = segList[iStU]->Yline_bbrrf();
00265 sx1 = segList[iStU]->Xline_slope();
00266 sy1 = segList[iStU]->Yline_slope();
00267 if (4 == m_debug) std::cout prt(xl1_0) prt(yl1_0) prt(sx1) prt(sy1) prt(omegag) << std::endl;
00268 if (fabs(omegag) >Constants::epsilon) {
00269 zintAU = findz_cyl(iAx, iStU, segList);
00270 xl1 = xl1_0 + sx1*zintAU;
00271 yl1 = yl1_0 + sy1*zintAU;
00272 rl1 = sqrt(xl1*xl1 + yl1*yl1);
00273 phiStU = atan2(yl1-yc, xl1-xc);
00274 if (4 == m_debug) cout prt(zintAU) prt(xl1) prt(yl1) prt(rl1) prt(phiStU) <<endl;
00275 } else {
00276 zintAU = -9999.;
00277 }
00278
00279
00280 zintAU_sl = findz_sl(iAx, iStU, segList);
00281 xl1_sl = xl1_0 + sx1*zintAU_sl;
00282 yl1_sl = yl1_0 + sy1*zintAU_sl;
00283 rl1_sl = sqrt(xl1_sl*xl1_sl + yl1_sl*yl1_sl);
00284 if (4 == m_debug) cout prt(zintAU_sl) prt(xl1_sl) prt(yl1_sl) prt(rl1_sl) <<endl;
00285
00286
00287 if ( (zintAU < -MdcxParameters::maxMdcZLen) &&
00288 (zintAU_sl < -MdcxParameters::maxMdcZLen) ) {
00289 if (4 == m_debug) std::cout << " SKIP by zintAU < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl;
00290 continue;
00291 }
00292 if ( (zintAU > MdcxParameters::maxMdcZLen) &&
00293 (zintAU_sl > MdcxParameters::maxMdcZLen) ) {
00294 if (4 == m_debug) std::cout << " SKIP by zintAU > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl;
00295 continue;
00296 }
00297
00298
00299
00300
00301 for (int iStV = -1; true; ) {
00302 if(!b_useGood_stV && (iStV >= (nSeg-1))) {
00303 if (4 == m_debug) std::cout << iStV <<" no segments in V slayer "<< std::endl;
00304 break;
00305 }
00306 if(b_useGood_stV && (iStV >= (nSeg-1))){
00307 iStV = -1;
00308 b_useGood_stV = false;
00309 }
00310 iStV++;
00311 if ((segList[iStV]->GetUsedOnHel()!=0)
00312 || (segList[iStU]->GetUsedOnHel()!=0)
00313 || (segList[iAx]->GetUsedOnHel()!=0)
00314 || (segList[iStV]->Fail() != 0)
00315 || (segList[iStV]->SuperLayer(0) != stVlay)
00316 || (segList[iStV]->Origin() != -1)
00317
00318 || (b_useGood_stV && (segList[iStV]->Nhits()< 4))
00319 ){ continue; }
00320
00321 if (4 == m_debug){
00322 std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl;
00323 std::cout<< "---------3.St V seg ------No."<<iStV<< std::endl;
00324 segList[iStV]->printSeg();
00325 std::cout<<" omega of seg iStV = "<<segList[iStV]->Omega()<<std::endl;
00326 }
00327
00328
00329 const HepAList<MdcxHit>& hitsStV = segList[iStV]->XHitList();
00330 double dPhiAV = fabs(hitsAx[0]->phiMid() - hitsStV[0]->phiMid());
00331 if (dPhiAV > Constants::pi) dPhiAV = Constants::twoPi - dPhiAV;
00332
00333 if(4 == m_debug) {
00334 if (dPhiAV > max_dPhiAV){
00335 cout << "SKIP by dPhiAV " <<dPhiAV <<" > "<< max_dPhiAV << endl;
00336 }else{
00337 cout << "KEEP by dPhiAV " <<dPhiAV<< " <= " << max_dPhiAV << endl;
00338 }
00339 std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
00340 <<" StV "<<hitsStV[0]->phiMid() << std::endl;
00341 std::cout<< std::endl;
00342 }
00343 if (g_dPhiAV) g_dPhiAV->fill(dPhiAV,segList[iStV]->SuperLayer());
00344 if (dPhiAV > max_dPhiAV) { continue; }
00345
00346
00347
00348 xl2_0 = segList[iStV]->Xline_bbrrf();
00349 yl2_0 = segList[iStV]->Yline_bbrrf();
00350 sx2 = segList[iStV]->Xline_slope();
00351 sy2 = segList[iStV]->Yline_slope();
00352 if (fabs(omegag) > Constants::epsilon) {
00353 zintAV = findz_cyl(iAx, iStV, segList);
00354 xl2 = xl2_0 + sx2*zintAV;
00355 yl2 = yl2_0 + sy2*zintAV;
00356 rl2 = sqrt(xl2*xl2 + yl2*yl2);
00357 if (4 == m_debug){
00358 segList[iAx]->printSeg();
00359 segList[iStV]->printSeg();
00360 cout << "zintAV " << zintAV << endl;
00361 }
00362 phiStV = atan2(yl2-yc, xl2-xc);
00363 } else {
00364 zintAV = -9999.;
00365 }
00366
00367
00368 zintAV_sl = findz_sl(iAx, iStV, segList);
00369 xl2_sl = xl2_0 + sx2*zintAV_sl;
00370 yl2_sl = yl2_0 + sy2*zintAV_sl;
00371 rl2_sl = sqrt(xl2_sl*xl2_sl + yl2_sl*yl2_sl);
00372
00373
00374 if ( (zintAV < -MdcxParameters::maxMdcZLen) &&
00375 (zintAV_sl < -MdcxParameters::maxMdcZLen) ) {
00376 if (4 == m_debug) std::cout << " SKIP by zintAV < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl;
00377 continue;
00378 }
00379 if ( (zintAV > MdcxParameters::maxMdcZLen) &&
00380 (zintAV_sl > MdcxParameters::maxMdcZLen) ) {
00381 if (4 == m_debug) std::cout << " SKIP by zintAV > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl;
00382 continue;
00383 }
00384
00385
00386 if (4 == m_debug) {
00387
00388 segList[iAx]->printSeg();
00389 segList[iStU]->printSeg();
00390 segList[iStV]->printSeg();
00391 }
00392 MdcxFittedHel fithel;
00393 double helixFitSigma = MdcxParameters::helixFitSigma;
00394 double first_prob = 0.0;
00395 HepAList<MdcxHit> hitlist;
00396 hitlist.append(hitsAx);
00397 hitlist.append(hitsStU);
00398 hitlist.append(hitsStV);
00399 if(g_trkllmk) for (int i=0; i<hitlist.length(); i++ ){ g_trkllmk->fill(hitlist[i]->Layer()); }
00400
00401
00402
00403
00404
00405 floop = 1;
00406 while (floop) {
00407 if(4 == m_debug) std::cout<< "---------4.Ax circle fit------"<< std::endl;
00408 if (fabs(omegag) < Constants::epsilon) {
00409 if(4 == m_debug) std::cout<<"SKIP by omegag==0 "<<std::endl;
00410 break;
00411 }
00412 if ( (zintAU < -MdcxParameters::maxMdcZLen) ||
00413 (zintAU > MdcxParameters::maxMdcZLen) ){
00414 if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
00415 break;
00416 }
00417 if ( (zintAV < -MdcxParameters::maxMdcZLen) ||
00418 (zintAV > MdcxParameters::maxMdcZLen) ){
00419 if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
00420 break;
00421 }
00422
00423
00424
00425 if(4 == m_debug) cout<< "dlen calc. slay U "<<segList[iStU]->SuperLayer()<< " slay V "<<segList[iStV]->SuperLayer()<<endl;
00426 double fltLenUV = dlen(segList[iStU]->SuperLayer(), phiStU,segList[iStV]->SuperLayer(), phiStV, omegag);
00427
00428 if(4 == m_debug){
00429 cout prt(fltLenUV) prt(phiStU) prt(phiStV) prt(omegag)<< std::endl;
00430 }
00431 if (fltLenUV > MdcxParameters::maxDlen) {
00432 if(4 == m_debug) std::cout<<"SKIP by fltLenUV > "<<MdcxParameters::maxDlen<<std::endl;
00433 break;
00434 }
00435 tanlg = (zintAV - zintAU)/fltLenUV;
00436 if(4 == m_debug) cout<< "dlen calc. slay A "<<segList[iAx]->SuperLayer()<< " slay U "<<segList[iStU]->SuperLayer()<<endl;
00437 double fltLenAU = dlen(segList[iAx]->SuperLayer(), phiAx, segList[iStU]->SuperLayer(), phiStU, omegag);
00438
00439 if(m_xtupleSegComb){
00440 m_segCombOmega = omegag;
00441 m_segCombSameAU = testFromSameTrack(segList[iAx],segList[iStU]);
00442 m_segCombSameUV = testFromSameTrack(segList[iStU],segList[iStV]);
00443 m_segCombSlA= segList[iAx]->SuperLayer();
00444 m_segCombSlU= segList[iStU]->SuperLayer();
00445 m_segCombSlV= segList[iStV]->SuperLayer();
00446 m_segCombPhiA = phiAx;
00447 m_segCombPhiU = phiStU;
00448 m_segCombPhiV = phiStV;
00449 m_segCombDLenUV = fltLenUV;
00450 m_segCombDLenAU = fltLenAU;
00451 m_xtupleSegComb->write();
00452 }
00453
00454 z0g = zintAU - fltLenAU*tanlg;
00455 if (4 == m_debug) cout prt(z0g) prt(tanlg) prt(fltLenAU) prt(zintAU) prt(zintAV)<< endl;
00456
00457
00458 if ((z0g < -MdcxParameters::maxMdcZLen) || (z0g > MdcxParameters::maxMdcZLen)) {
00459 if (4 == m_debug) std::cout<<"SKIP by z0g out of range "<<z0g<<">"<<MdcxParameters::maxMdcZLen<<std::endl;
00460 break;
00461 }
00462 ztest = z0g + fltAx*tanlg;
00463 if (4 == m_debug) std::cout prt(ztest) prt(fltAx)<<std::endl;
00464
00465
00466
00467 MdcxHel ghel(segList[iAx]->D0(), segList[iAx]->Phi0(), segList[iAx]->Omega(),
00468 ztest, tanlg, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref());
00469
00470 MdcxFittedHel firstfit(hitlist,ghel,helixFitSigma);
00471 first_prob = firstfit.Prob();
00472 if (4 == m_debug) {
00473 std::cout<<" after first fit: "<<std::endl;
00474 firstfit.FitPrint();
00475 }
00476 if (first_prob >= probmin) {
00477 MdcxHel helixOrigin= TakeToOrigin(firstfit);
00478 MdcxFittedHel nextfit(hitlist,helixOrigin,helixFitSigma);
00479 first_prob = nextfit.Prob();
00480 fithel=nextfit;
00481 if (g_trklgood) for (int i=0; i<nextfit.Nhits(); i++){ g_trklgood->fill(nextfit.Layer(i));}
00482 }
00483 if ( g_trklfirstProb) g_trklfirstProb->fill(first_prob);
00484 if (4 == m_debug) cout << " prob after helix fitting = " << first_prob << endl;
00485 floop = 0;
00486 }
00487
00488
00489
00490
00491
00492 floop = 1;
00493 while (floop) {
00494 if (4 == m_debug) std::cout<< "---------4.2 straight line fit------"<< std::endl;
00495 if (4 == m_debug) cout << " zintAU_sl " << zintAU_sl << " zintAV_sl " << zintAV_sl << endl;
00496 if (first_prob > probmin) {
00497 if (4 == m_debug) cout << " helix fit good , exit straight line fit." << endl;
00498 break;
00499 }
00500
00501 if ( (zintAU_sl < -MdcxParameters::maxMdcZLen) ||
00502 (zintAU_sl > MdcxParameters::maxMdcZLen) ) break;
00503 if ( (zintAV_sl < -MdcxParameters::maxMdcZLen) ||
00504 (zintAV_sl > MdcxParameters::maxMdcZLen) ) break;
00505
00506 double dx = xl2_sl - xl1_sl;
00507 double dy = yl2_sl - yl1_sl;
00508 double dl = sqrt(dx*dx + dy*dy);
00509 tanlg_sl = (zintAV_sl - zintAU_sl)/dl;
00510 dx = xl1_sl + d0g_sl*sin(phi0g_sl);
00511 dy = yl1_sl - d0g_sl*cos(phi0g_sl);
00512 dl = sqrt(dx*dx + dy*dy);
00513 z0g_sl = zintAU_sl - dl*tanlg_sl;
00514 if (4 == m_debug) cout prt(d0g_sl) prt(phi0g_sl) prt(z0g_sl) prt(tanlg_sl)<< endl;
00515
00516 if ((z0g_sl < -MdcxParameters::maxMdcZLen) || (z0g_sl > MdcxParameters::maxMdcZLen)){
00517 if(4 == m_debug) std::cout << "SKIP by z0g_sl:"<<z0g_sl<<" > "<<MdcxParameters::maxMdcZLen << std::endl;
00518 break;
00519 }
00520
00521 ztest_sl = z0g_sl + fltL_sl*tanlg_sl;
00522 if (4 == m_debug) cout prt(ztest_sl) << endl;
00523
00524
00525 MdcxHel ghel(segList[iAx]->D0_sl_approx(), segList[iAx]->Phi0_sl_approx(), omegag_sl,
00526 ztest_sl, tanlg_sl, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref());
00527
00528 MdcxFittedHel firstfit(hitlist, ghel, helixFitSigma);
00529 first_prob = firstfit.Prob();
00530
00531 if(4 == m_debug){
00532 ghel.print();
00533 std::cout <<"first_prob "<<first_prob << std::endl;
00534 }
00535
00536 if (first_prob >= probmin) {
00537 MdcxHel helixOrigin = TakeToOrigin(firstfit);
00538 MdcxFittedHel nextfit(hitlist, helixOrigin, helixFitSigma);
00539 first_prob = nextfit.Prob();
00540 fithel = nextfit;
00541 if (4 == m_debug) {
00542 cout << " three seg sl nextfit" << endl;
00543 nextfit.FitPrint();
00544 }
00545 }
00546 floop = 0;
00547 }
00548
00549
00550
00551
00552 floop = 1;
00553 while (floop) {
00554 floop = 0;
00555 if(4 == m_debug)std::cout<< "---------5. try add seg to helix------"<< std::endl;
00556 if (first_prob < probmin) {
00557 if(4 == m_debug) std::cout<< "prob"<<first_prob<< " "
00558 <<probmin<<" fit NOT good, exit add segs "<< std::endl;
00559 break;
00560 }
00561 float bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
00562 int bndof = fithel.Nhits() - 5;
00563 float bprob = Mdcxprobab(bndof, bchisq);
00564 trkCount++;
00565 segList[iAx]->SetUsedOnHel(trkCount);
00566 segList[iStU]->SetUsedOnHel(trkCount);
00567 segList[iStV]->SetUsedOnHel(trkCount);
00568
00569 if (4 == m_debug){
00570 cout << "in add seg to trail helix, klm seg:" << endl;
00571 segList[iAx]->printSeg();
00572 segList[iStU]->printSeg();
00573 segList[iStV]->printSeg();
00574 }
00575
00576
00577 for (int iSlay = 0; iSlay < Constants::nSuperLayer; iSlay++) {
00578 int ilay = MdcxParameters::layerSet2AddSeg[iSegCombo][iSlay];
00579 for (int iSeg = 0; iSeg < nSeg; iSeg++) {
00580 if ((segList[iSeg]->SuperLayer(0) != ilay)
00581 ||(segList[iSeg]->GetUsedOnHel() != 0)
00582 || (segList[iSeg]->Fail() != 0)
00583 || (segList[iSeg]->Origin() != -1)) continue;
00584 if(4 == m_debug) {
00585 std::cout<< endl<< "== add seg "<< iSeg<<" ";
00586 segList[iSeg]->printSeg();
00587 }
00588
00589
00590 float phiAxSeg = segList[iAx]->XHitList()[0]->phiMid();
00591 float phiAddSeg = segList[iSeg]->XHitList()[0]->phiMid();
00592 double phiKNSegDiff = fabs(phiAxSeg - phiAddSeg);
00593
00594
00595
00596
00597 int t_same=-999;
00598 if(m_xtupleAddSeg1){
00599
00600 t_same = testFromSameTrack(segList[iAx],segList[iSeg]);
00601 m_addSegSame = t_same;
00602 m_addSegSeedSl = segList[iAx]->SuperLayer();
00603 m_addSegSeedPhi = segList[iAx]->XHitList()[0]->phiMid();
00604 m_addSegSeedPhiLay = segList[iAx]->XHitList()[0]->Layer();
00605 m_addSegSeedPhi0 = segList[iAx]->Phi0_sl_approx();
00606 m_addSegSeedD0 = segList[iAx]->D0_sl_approx();
00607 m_addSegAddSl = segList[iSeg]->SuperLayer();
00608 m_addSegAddPhi = segList[iSeg]->XHitList()[0]->phiMid();
00609 m_addSegAddPhiLay = segList[iSeg]->XHitList()[0]->Layer();
00610 m_addSegAddPhi0 = segList[iSeg]->Phi0_sl_approx();
00611 m_addSegAddD0 = segList[iSeg]->D0_sl_approx();
00612 m_xtupleAddSeg1->write();
00613 }
00614
00615 if (phiKNSegDiff > 0.5*Constants::pi && phiKNSegDiff < 1.5*Constants::pi )
00616
00617 {
00618 if(4 == m_debug) std::cout<< " SKIP by phi diff,same "<< t_same <<" Ax "<<
00619 phiAxSeg<<" iSeg "<<phiAddSeg<<" diff "<<phiKNSegDiff << std::endl;
00620 continue;
00621 }else{
00622 if(4 == m_debug) std::cout<< " phi diff OK, Ax "<<
00623 phiAxSeg<<" added "<<phiAddSeg<<" diff="<<phiKNSegDiff << std::endl;
00624 }
00625
00626
00627 xp = segList[iSeg]->Xref() - segList[iSeg]->SinPhi0()*segList[iSeg]->D0();
00628 yp = segList[iSeg]->Yref() + segList[iSeg]->CosPhi0()*segList[iSeg]->D0();
00629 const HepAList<MdcxHit>& hitsSegAdd= segList[iSeg]->XHitList();
00630 double proca = fithel.Doca( hitsSegAdd[0]->wx(), hitsSegAdd[0]->wy(),
00631 hitsSegAdd[0]->wz(), xp, yp );
00632 if (g_trklproca) g_trklproca->fill(proca);
00633
00634
00635 if(m_xtupleAddSeg2){
00636 m_addSegEvtNo = g_eventNo;
00637 m_addSegPoca = proca;
00638 m_addSegSlayer = iSlay;
00639 m_addSegLen = fithel.Doca_Len();
00640 m_xtupleAddSeg2->write();
00641 }
00642
00643 if (!((fabs(proca) < MdcxParameters::maxProca)&&(fithel.Doca_Len()<fithel.Lmax())) ){
00644
00645 if(4 == m_debug){
00646 std::cout<< " SKIP by fabs(proca):"<< fabs(proca)<< "<"<<MdcxParameters::maxProca <<" or Doca_Len:"<<fithel.Doca_Len() <<"<"<<fithel.Lmax()<< std::endl;
00647 }
00648 }else{
00649 if(4 == m_debug){
00650 std::cout<< " proca & len OK, |proca| "<< fabs(proca)<< "<"<<MdcxParameters::maxProca <<" && Doca_Len:"<<fithel.Doca_Len() <<"<"<<fithel.Lmax()<< std::endl;
00651 }
00652 }
00653 if(fithel.Doca_Len()<0) continue;
00654
00655 if ((fabs(proca)<MdcxParameters::maxProca)&&(fithel.Doca_Len()<fithel.Lmax())) {
00656 MdcxFittedHel newhel;
00657 newhel.Grow(fithel, hitsSegAdd);
00658 if (newhel.Nhits() != hitlist.length()) {
00659 if (4 == m_debug){
00660 cout<<" rcs "<<newhel.Rcs()<<"<=?"<<MdcxParameters::maxRcsInAddSeg<< endl;
00661 }
00662
00663
00664 if (newhel.Prob() < probmin){
00665 if(4 == m_debug){
00666 cout << " prob " << newhel.Prob() << "<"<<probmin<<", ReFit" << endl;
00667 }
00668 newhel.ReFit();
00669 }
00670 if(4 == m_debug) {
00671 cout<<" refit prob "<<newhel.Prob()<<"<"<<probmin<<" Fail? "<<newhel.Fail()<< endl;
00672 }
00673
00674 if ( (newhel.Prob() >= probmin) && (newhel.Fail() == 0) ) {
00675
00676 fithel = newhel;
00677 hitlist = newhel.XHitList();
00678 segList[iSeg]->SetUsedOnHel(trkCount);
00679 bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
00680 bndof = fithel.Nhits() - 5;
00681 float tprob = Mdcxprobab(bndof, bchisq);
00682 if (tprob > bprob) bprob = tprob;
00683 if (4 == m_debug) {
00684 cout << " segment ADDED, prob=" << newhel.Prob() << endl;
00685 }
00686 }
00687 else{
00688 if(4 == m_debug)std::cout<<" fit FAILED"<<std::endl;
00689 }
00690 } else {
00691 segList[iSeg]->SetUsedOnHel(trkCount);
00692 if (4 == m_debug) cout << " segment ADDED, but no new hits" << endl;
00693 }
00694 }
00695 }
00696 }
00697 if (((fithel.Prob() < 0.90) && (fithel.Nhits() == 12)) || (fithel.Fail() != 0)) {
00698 for (int i = 0; i < nSeg; i++) { if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0); }
00699 trkCount--;
00700 break;
00701 }
00702
00703 if(4 == m_debug){
00704 std::cout<< "Track after add segs" << std::endl;
00705 fithel.FitPrint();
00706 }
00707
00708
00709 if(4 == m_debug) std::cout<< "---------6. flip------"<< std::endl;
00710 fithel.flip();
00711
00712
00713 int kki = 0;
00714 int notduplicate = 1;
00715 while(MdcxTrklist[kki]) {
00716 if (4 == m_debug)std::cout<< "---------7. test saved track------"<< std::endl;
00717 const HepAList<MdcxHit>& junk=MdcxTrklist[kki]->XHitList();
00718 int overl = 0;
00719 int kkj = 0;
00720 while (junk[kkj]) {
00721 int kkl = 0;
00722 while (hitlist[kkl]) {
00723 if (hitlist[kkl] == junk[kkj]) overl++;
00724 kkl++;
00725 }
00726 kkj++;
00727 }
00728 if (4 == m_debug) cout << "overlap with track # " << kki << " is "
00729 << junk.length() << " hitList " << hitlist.length() << " overlap " << overl << endl;
00730 if ( (hitlist.length()/4.) <= overl) notduplicate = 0;
00731
00732
00733 kki++;
00734 }
00735
00736 if(g_trklhelix) for (int iHit=0; iHit<fithel.Nhits(); iHit++){ g_trklhelix->fill(fithel.Layer(iHit)); }
00737
00738
00739 if (notduplicate) {
00740 if(4 == m_debug) std::cout<< "---------8. test worth saving?------"<< std::endl;
00741 MdcxFittedHel *finehel = new MdcxFittedHel(hitlist, fithel);
00742
00743 if(4 == m_debug) {
00744 std::cout<<__FILE__
00745 <<" finehel Prob>0.001 "<<finehel->Prob()
00746 <<" nhits "<<finehel->Nhits() <<">=15? "
00747 <<" bprob "<< bprob
00748 <<" >=?probmin "<< probmin
00749 <<" Fail==0? "<< finehel->Fail()
00750 << std::endl;
00751 finehel->FitPrint();
00752 }
00753
00754 if (( (finehel->Prob()>0.001) || (finehel->Nhits()>=15) || (bprob>probmin)
00755
00756 ) && (finehel->Fail() == 0) ) {
00757 nfound++;
00758 sumprob += finehel->Prob();
00759
00760
00761
00762
00763 int nodrop = 0;
00764 if ((finehel->Prob() > 0.05) || nodrop ) {
00765 MdcxTrklist.append(finehel);
00766
00767 for (int iHit=0; iHit<finehel->Nhits(); iHit++){
00768 if(g_trklappend1) g_trklappend1->fill(finehel->Layer(iHit));
00769 }
00770
00771 } else {
00772 MdcxFittedHel* drophel = new MdcxFittedHel( drophits(finehel) );
00773 float ratdrop = float(drophel->Nhits()) / float(finehel->Nhits());
00774 if (4 == m_debug) cout << "APPEND track " << trkCount << ", drops "
00775 << finehel->Nhits()-drophel->Nhits()
00776 << " helixnhit "<<finehel->Nhits()<< " drophit "<<drophel->Nhits()
00777 << " hits, drop rate="<<ratdrop <<">?"<<0.5
00778 << " prob "<<drophel->Prob() <<" >?"<<finehel->Prob()
00779 << " fail==0? "<<drophel->Fail()
00780 << endl;
00781 if ((drophel->Prob() > finehel->Prob()) &&
00782 (ratdrop > 0.50) && (drophel->Fail() == 0)) {
00783 if (4 == m_debug) cout << "append drop helix " << endl;
00784 MdcxTrklist.append(drophel);
00785 if(g_trklappend2) for (int iHit=0;iHit<drophel->Nhits();iHit++){g_trklappend2->fill(drophel->Layer(iHit));}
00786 delete finehel;
00787 } else {
00788 if (4 == m_debug) cout << "append fine helix " << endl;
00789 MdcxTrklist.append(finehel);
00790 if(g_trklappend3)for(int iHit=0;iHit<finehel->Nhits();iHit++){g_trklappend3->fill(finehel->Layer(iHit));}
00791 delete drophel;
00792 }
00793 }
00794
00795
00796 int nwl = MdcxTrklist.length() - 1;
00797 if ((4 == m_debug) && (nwl > -1)) {
00798 if (4 == m_debug) {
00799 cout << endl << "found track " << trkCount<<std::endl;
00800 MdcxTrklist[nwl]->print();
00801 MdcxTrklist[nwl]->FitPrint();
00802 }
00803 }
00804 } else {
00805 for (int i = 0; i < nSeg; i++) {
00806 if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0);
00807 }
00808 delete finehel;
00809 trkCount--;
00810 }
00811 }
00812 }
00813 }
00814 }
00815 }
00816 if(4==m_debug) cout<< "end of this combo"<<endl;
00817 }
00818 if(4== m_debug) cout << " In MdcxFindTracks, found " << trkCount << " good tracks" << endl;
00819 }
00820
00821 void MdcxFindTracks::print(ostream &o) {
00822 o<<"ngood:"<<ngood<<endl;
00823 o<<"npure:"<<npure<<endl;
00824 o<<"nbad:"<<nbad<<endl;
00825 o<<"nfail:"<<nfail<<endl;
00826 o<<"nhitsmc:"<<nhitsmc<<endl;
00827 o<<"nhitsgd:"<<nhitsgd<<endl;
00828 }
00829
00830 void MdcxFindTracks::plot()const {
00831 }
00832
00833 void MdcxFindTracks::printpar(ostream &o) {
00834 o << "|MdcxFindTracks| parameters" << endl;
00835 o << "probmin:" << probmin << endl;
00836 o << "curvmax:" << curvmax << endl;
00837 }
00838
00839 void MdcxFindTracks::beginmore(){
00840 ngood=0; nbad=0; nfail=0; npure=0;
00841 nhitsmc=0; nhitsgd=0;
00842 printpar(cout);
00843 }
00844
00845 void MdcxFindTracks::endmore(){
00846 print(cout);
00847 }
00848
00849 double MdcxFindTracks::findz_sl(int iAx, int iStU, const HepAList<MdcxSeg> &segList){
00850 double Ap = -sin(segList[iAx]->Phi0_sl_approx());
00851 double Bp = cos(segList[iAx]->Phi0_sl_approx());
00852 double xp = segList[iAx]->Xref() + Ap*segList[iAx]->D0_sl_approx();
00853 double yp = segList[iAx]->Yref() + Bp*segList[iAx]->D0_sl_approx();
00854 double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0();
00855 double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
00856
00857 const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
00858 double Cl = hitsStU[0]->wz();
00859 double ndotl = Bp*hitsStU[0]->wy() + Ap*hitsStU[0]->wx();
00860
00861
00862
00863
00864 double zint = (ndotl != 0.0) ? ((Bp*(yp-yl)+Ap*(xp-xl))*Cl/ndotl) : -1000.;
00865
00866 return zint;
00867 }
00868
00869 double MdcxFindTracks::findz_cyl(int iAx, int iStU, const HepAList<MdcxSeg> &segList) {
00870 double zint = -1000.0;
00871 double r = 1.0/fabs(segList[iAx]->Omega());
00872 double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0();
00873 double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
00874 double sx = segList[iStU]->Xline_slope();
00875 double sy = segList[iStU]->Yline_slope();
00876 double cx = xl_0 - segList[iAx]->Xc();
00877 double cy = yl_0 - segList[iAx]->Yc();
00878 double a = sx*sx + sy*sy;
00879 double b = 2.0*(sx*cx + sy*cy);
00880
00881 double bsq = b * b;
00882 double fac = 4.0 * a * (cx*cx + cy*cy - r*r);
00883 double ta = 2.0 * a;
00884 if (fac < bsq) {
00885 double sfac = sqrt(bsq - fac);
00886 double z1 = -(b - sfac)/ta;
00887 double z2 = -(b + sfac)/ta;
00888 zint = (fabs(z1) < fabs(z2)) ? z1 : z2;
00889 }
00890
00891 return zint;
00892 }
00893
00894 double MdcxFindTracks::dlen(int slay1, double p1, int slay2, double p2, double om){
00895
00896 double dp;
00897 if((slay1==2 || slay1==3 || slay1 == 4 || slay1 == 9 || slay1 ==10) ){
00898 dp = p2 - p1;
00899 }else if((slay2 >= slay1)) {
00900 dp = p2 - p1;
00901 }else {
00902 dp = p1 - p2;
00903 }
00904
00905 if (om < 0.0) {
00906 while (dp < -Constants::pi) dp += Constants::pi;
00907 while (dp > 0.0) dp -= 1*Constants::pi;
00908 } else {
00909 while (dp < 0.0) dp += 1*Constants::pi;
00910 while (dp > Constants::pi) dp -= 1*Constants::pi;
00911 }
00912
00913 double dl = dp/om;
00914
00915
00916
00917
00918 return dl;
00919 }
00920
00921 MdcxFittedHel MdcxFindTracks::drophits(MdcxFittedHel* finehel) {
00922 MdcxHel mdcxHel = (MdcxHel) *finehel;
00923 const HepAList<MdcxHit>& hitlist = finehel->XHitList();
00924 HepAList<MdcxHit> nwhitlist = hitlist;
00925 int nremove = 0;
00926 int kkk = 0;
00927 while(nwhitlist[kkk]) {
00928 double pull = nwhitlist[kkk]->pull(mdcxHel);
00929 int layer = nwhitlist[kkk]->Layer();
00930
00931 if(m_xtupleDropHits){
00932 float t_doca = mdcxHel.Doca(*(nwhitlist[kkk]));
00933 float t_e = nwhitlist[kkk]->e(t_doca);
00934 m_segDropHitsEvtNo = g_eventNo;
00935 m_segDropHitsLayer = layer;
00936 m_segDropHitsWire = nwhitlist[kkk]->WireNo();
00937 m_segDropHitsPull = pull;
00938 m_segDropHitsDoca = t_doca;
00939 m_segDropHitsSigma = t_e;
00940 m_segDropHitsDrift = nwhitlist[kkk]->d(mdcxHel);
00941 m_segDropHitsMcTkId = nwhitlist[kkk]->getDigi()->getTrackIndex();
00942 m_xtupleDropHits->write();
00943 }
00944
00945 if (fabs(pull) > MdcxParameters::dropHitsSigma[layer]){
00946 nwhitlist.remove(kkk);
00947 nremove++;
00948 } else {
00949 kkk++;
00950 }
00951 }
00952 if(m_debug==4) cout<< " MdcxFindTracks::drophits drop "<<nremove<<" hits "
00953 << " nhits of drop helix = "<<nwhitlist.length()<<endl;
00954
00955 MdcxFittedHel newhel(nwhitlist,mdcxHel);
00956 return newhel;
00957 }
00958
00959 void MdcxFindTracks::resout(MdcxFittedHel* finehel) {
00960 MdcxHel mdcxHel = (MdcxHel) *finehel;
00961 HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
00962 int kkk = 0; while(hitlist[kkk]){hitlist[kkk]->SetConstErr(0); kkk++;}
00963 MdcxFittedHel newhel(hitlist, mdcxHel);
00964 if(4 == m_debug) newhel.FitPrint(mhel,cout);
00965 kkk = 0;
00966 while(hitlist[kkk]) {
00967 float doca = newhel.Doca(*hitlist[kkk]);
00968 float zh = newhel.Doca_Zh();
00969 float tof = newhel.Doca_Tof();
00970 float dd = hitlist[kkk]->d(newhel);
00971 float tt = hitlist[kkk]->tcor(zh, tof);
00972 if (4 == m_debug){
00973 cout << "MdcxFindTracks::resout ("
00974 << hitlist[kkk]->Layer() << ","
00975 << hitlist[kkk]->WireNo() << ") dt "
00976 << tt << " resid "
00977 << fabs(dd)-fabs(doca) << " doca "
00978 << fabs(doca) << endl;
00979 }
00980 kkk++;
00981 }
00982 kkk=0;
00983 while(hitlist[kkk]){hitlist[kkk]->SetConstErr(1); kkk++;}
00984 }
00985
00986 MdcxHel MdcxFindTracks::TakeToOrigin(MdcxHel& ihel) {
00987 double lng;
00988 double omegag = ihel.Omega();
00989 double phi0g = ihel.Phi0();
00990 double phi0s = phi0g;
00991 double xp = ihel.Xref() - ihel.SinPhi0()*ihel.D0();
00992 double yp = ihel.Yref() + ihel.CosPhi0()*ihel.D0();
00993 double d0g = yp*cos(phi0g) - xp*sin(phi0g);
00994 if (omegag != 0.0) {
00995 phi0g = phi0g - Constants::pi;
00996 double xc = sin(phi0g)/omegag;
00997 double yc = -cos(phi0g)/omegag;
00998 double t1 = -xc - xp;
00999 double t2 = yc + yp;
01000 double phinew = atan2(t1, t2);
01001 if (omegag > 0.0) phinew += Constants::pi;
01002 if (phinew > Constants::pi) phinew -= Constants::twoPi;
01003 double xh = xc - sin(phinew)/omegag;
01004 double yh = yc + cos(phinew)/omegag;
01005 double x0 = xh + xp;
01006 double y0 = yh + yp;
01007 d0g = sqrt(x0*x0 + y0*y0);
01008 phi0g = phinew + Constants::pi;
01009 double sd0 = x0*sin(phi0g) - y0*cos(phi0g);
01010 if (sd0 > 0.0) d0g = -d0g;
01011 lng = dlen(-2, phi0g, -1, phi0s, omegag);
01012
01013 } else {
01014 lng = ihel.Xref()*ihel.CosPhi0() + ihel.Yref()*ihel.SinPhi0();
01015 }
01016 double tanlg = ihel.Tanl();
01017 double z0g = ihel.Z0() - lng*tanlg;
01018
01019 double t0g = ihel.T0();
01020 int codeg = ihel.Code();
01021 MdcxHel ohel(d0g, phi0g, omegag, z0g, tanlg, t0g, codeg);
01022 return ohel;
01023 }
01024
01025 bool MdcxFindTracks::testFromSameTrack(MdcxSeg* seg1, MdcxSeg* seg2){
01026
01027
01028
01029
01030
01031
01032
01033
01034
01035
01036
01037
01038
01039 int trackId = -9999;
01040 for (int i =0; i<seg1->Nhits(); i++){
01041 int tkid = seg1->XHitList()[i]->getDigi()->getTrackIndex();
01042 if(tkid>=1000) tkid -=1000;
01043 if((trackId == -9999)){
01044 if(tkid>=0) trackId = tkid;
01045 }else{
01046 if(tkid!=trackId) return false;
01047 }
01048 }
01049 for (int i =0; i<seg2->Nhits(); i++){
01050 int tkid = seg2->XHitList()[i]->getDigi()->getTrackIndex();
01051 if(tkid>=1000) tkid -=1000;
01052 if((trackId == -9999)){
01053 return false;
01054 }else{
01055 if(tkid!=trackId) return false;
01056 }
01057 }
01058 return true;
01059 }