Main Page | Namespace List | Class Hierarchy | Alphabetical List | Class List | Directories | File List | Namespace Members | Class Members | File Members | Related Pages

MdcxFindTracks Class Reference

#include <MdcxFindTracks.h>

List of all members.

Public Member Functions

const HepAList< MdcxFittedHel > & GetMdcxTrklist ()
const HepAList< MdcxFittedHel > & GetMdcxTrklist ()
 MdcxFindTracks (const HepAList< MdcxSeg > &f, int debug=0)
 MdcxFindTracks ()
 MdcxFindTracks (const HepAList< MdcxSeg > &f, int debug=0)
 MdcxFindTracks ()
void plot () const
void plot () const
void print (std::ostream &o)
void print (std::ostream &o)
void process (const HepAList< MdcxSeg > &f)
void process (const HepAList< MdcxSeg > &f)
virtual ~MdcxFindTracks ()
virtual ~MdcxFindTracks ()

Protected Member Functions

void beginmore ()
void beginmore ()
double dlen (double p1, double p2, double om)
double dlen (double p1, double p2, double om)
MdcxFittedHel drophits (MdcxFittedHel *f)
MdcxFittedHel drophits (MdcxFittedHel *f)
void endmore ()
void endmore ()
double findz_cyl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
double findz_cyl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
double findz_sl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
double findz_sl (int i1, int i2, const HepAList< MdcxSeg > &slclus)
void KillList ()
void KillList ()
void printpar (std::ostream &o)
void printpar (std::ostream &o)
void resout (MdcxFittedHel *f)
void resout (MdcxFittedHel *f)
MdcxHel TakeToOrigin (MdcxHel &)
MdcxHel TakeToOrigin (MdcxHel &)

Protected Attributes

float curvmax
int m_debug
bool m_doSag
HepAList< MdcxFittedHel > MdcxTrklist
HepAList< MdcxFittedHel > MdcxTrklist
MdcxHel mhel
int nbad
int nfail
int ngood
int nhitsgd
int nhitsmc
int npure
float probmin

Constructor & Destructor Documentation

MdcxFindTracks::MdcxFindTracks (   ) 
  

00064                               { 
00065   probmin = MdcxParameters::minTrkProb;
00066   curvmax = 1000000.0;
00067 }

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

00069                                                                          { 
00070   probmin = MdcxParameters::minTrkProb;
00071   curvmax = 1000000.0;
00072   m_debug = debug;
00073   process(segList); 
00074 }

MdcxFindTracks::~MdcxFindTracks (   )  [virtual]
  

00076 { KillList(); }

MdcxFindTracks::MdcxFindTracks (   ) 
 

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

virtual MdcxFindTracks::~MdcxFindTracks (   )  [virtual]
 

Member Function Documentation

void MdcxFindTracks::beginmore (   )  [protected]
 

void MdcxFindTracks::beginmore (   )  [protected]
  

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

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

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

00784                                                           {
00785   double dp = p2 - p1;
00786   if (om < 0.0) {
00787     while (dp < -Constants::pi) dp += Constants::pi;
00788     while (dp > 0.0)   dp -= Constants::pi;
00789   } else {
00790     while (dp < 0.0)   dp += Constants::pi;
00791     while (dp > Constants::pi)  dp -= Constants::pi;
00792   }
00793   //double r = 1./om;
00794   //double dl = dp * r;
00795   double dl = dp/om;
00796   //if(4 == m_debug) cout prt(p1) prt(p2) prt(dp) prt(om) prt(r) prt(dl) << endl;
00797   return dl;
00798 }

MdcxFittedHel MdcxFindTracks::drophits (  MdcxFittedHel *  f  )  [protected]
 

MdcxFittedHel MdcxFindTracks::drophits (  MdcxFittedHel *  f  )  [protected]
  

00800                                                              {
00801   MdcxHel mdcxHel = (MdcxHel) *finehel;
00802   const HepAList<MdcxHit>& hitlist = finehel->XHitList();
00803   HepAList<MdcxHit> nwhitlist = hitlist;
00804   int kkk = 0;
00805   while(nwhitlist[kkk]) {
00806     double pull = nwhitlist[kkk]->pull(mdcxHel);
00807     float t_dd = nwhitlist[kkk]->d(mdcxHel);
00808     float t_doca = mdcxHel.Doca(*(nwhitlist[kkk]));
00809     float t_e =  nwhitlist[kkk]->e(t_doca);
00810     int layer = nwhitlist[kkk]->Layer();
00811     if(g_dropHitsSigma)g_dropHitsSigma->fill(layer,fabs(pull));
00812     if (fabs(pull) > MdcxParameters::dropHitsSigma[layer]){
00813       if (g_trkld) g_trkld->fill(t_dd);
00814       if (g_trkldl) g_trkldl->fill(t_dd,layer);
00815       if (g_trklel) g_trklel->fill(t_e,layer);
00816       if (g_trkld) g_trkld->fill(t_dd);
00817       if (g_trkldoca) g_trkldoca->fill(t_doca);
00818       if (g_trkllayer) g_trkldoca->fill(layer);
00819       nwhitlist.remove(kkk);
00820     } else {
00821       kkk++;
00822     }
00823   }
00824   MdcxFittedHel newhel(nwhitlist,mdcxHel);
00825   return newhel;
00826 }//endof drophits

void MdcxFindTracks::endmore (   )  [protected]
 

void MdcxFindTracks::endmore (   )  [protected]
  

00735                             {
00736   print(cout);
00737 }//endof endmore

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

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

00759                                                                                     {
00760   double zint = -1000.0;
00761   double r = 1.0/fabs(segList[iAx]->Omega());
00762   double xl_0 = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0(); 
00763   double yl_0 = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
00764   double sx = segList[iStU]->Xline_slope();
00765   double sy = segList[iStU]->Yline_slope();
00766   double cx = xl_0 - segList[iAx]->Xc();
00767   double cy = yl_0 - segList[iAx]->Yc();
00768   double a = sx*sx + sy*sy;
00769   double b = 2.0*(sx*cx + sy*cy);
00770   //double c = (cx*cx + cy*cy - r*r);
00771   double bsq = b * b;
00772   double fac = 4.0 * a * (cx*cx + cy*cy - r*r);
00773   double ta = 2.0 * a;
00774   if (fac < bsq) {
00775     double sfac = sqrt(bsq - fac);
00776     double z1   = -(b - sfac)/ta;
00777     double z2   = -(b + sfac)/ta;
00778     zint = (fabs(z1) < fabs(z2)) ? z1 : z2;
00779   }
00780   // cout << " findz_cyl segs " << iAx << " " << iStU << " zint " << zint << endl;
00781   return zint;
00782 }

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

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

00739                                                                                   {
00740   double Ap = -sin(segList[iAx]->Phi0_sl_approx());
00741   double Bp =  cos(segList[iAx]->Phi0_sl_approx());
00742   double xp = segList[iAx]->Xref() + Ap*segList[iAx]->D0_sl_approx();
00743   double yp = segList[iAx]->Yref() + Bp*segList[iAx]->D0_sl_approx();
00744   double xl = segList[iStU]->Xref() - segList[iStU]->SinPhi0()*segList[iStU]->D0(); 
00745   double yl = segList[iStU]->Yref() + segList[iStU]->CosPhi0()*segList[iStU]->D0();
00746 
00747   const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
00748   double Cl = hitsStU[0]->wz();
00749   double ndotl = Bp*hitsStU[0]->wy() + Ap*hitsStU[0]->wx();
00750   /*
00751      std::cout<<"ndotl   "<<ndotl<<" "<<Bp<<" "<<hitsStU[0]->wy()<<" "<<
00752      Ap<<" "<<hitsStU[0]->wx()<<std::endl;
00753      */
00754   double zint = (ndotl != 0.0) ? ((Bp*(yp-yl)+Ap*(xp-xl))*Cl/ndotl) : -1000.;
00755   // cout << " findz_sl segs " << iAx << " " << iStU << " zint " << zint << endl;
00756   return zint;
00757 }

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

00039 {return MdcxTrklist;}

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

00039 {return MdcxTrklist;}

void MdcxFindTracks::KillList (   )  [inline, protected]
  

00061 {HepAListDeleteAll(MdcxTrklist);}

void MdcxFindTracks::KillList (   )  [inline, protected]
  

00061 {HepAListDeleteAll(MdcxTrklist);}

void MdcxFindTracks::plot (   )  const
 

void MdcxFindTracks::plot (   )  const
  

00720                                {
00721 } // endof plot

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

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

00711                                      {
00712   o<<"ngood:"<<ngood<<endl;
00713   o<<"npure:"<<npure<<endl;
00714   o<<"nbad:"<<nbad<<endl;
00715   o<<"nfail:"<<nfail<<endl;
00716   o<<"nhitsmc:"<<nhitsmc<<endl;
00717   o<<"nhitsgd:"<<nhitsgd<<endl;
00718 } // endof print

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

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

00723                                         {
00724   o << "|MdcxFindTracks| parameters" << endl;
00725   o << "probmin:" << probmin << endl;
00726   o << "curvmax:" << curvmax << endl;
00727 }//endof printpar

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

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

ZOUJH

ZOUJH

ZOUJH

ZOUJH

ZOUJH

ZOUJH 

00078                                                              {
00079   static int   nfound  = 0;
00080   static float sumprob = 0.0; 
00081 
00082   MdcxHel thel(0., 0., 0., 0., 0.);
00083   mhel = thel;
00084   int nSeg = segList.length();
00085   for (int i = 0; i < nSeg ; i++) { segList[i]->SetUsedOnHel(0); }
00086 
00087   double xp, yp, xl1, yl1, rl1, xl2, yl2, rl2;
00088   double d0g, phi0g, omegag, z0g, tanlg;
00089   double d0g_sl, phi0g_sl, omegag_sl, z0g_sl, tanlg_sl;
00090   double zintAU, zintAV, zintAU_sl, zintAV_sl = 9999.;
00091   double rl1_sl, rl2_sl;
00092   double xc, yc;
00093   double sinPhi0g_sl, cosPhi0g_sl, xp_sl, yp_sl;
00094   double phiAx, phiStU, phiStV, xl1_sl, yl1_sl, xl2_sl, yl2_sl, xl1_0, yl1_0, xl2_0, yl2_0;
00095   double sx1, sy1, sx2, sy2;
00096   double x0g, y0g;  
00097   double fltAx, ztest, fltL_sl, ztest_sl;  
00098   int floop, trkCount= 0;
00099 
00100   //The main loop; first get the 3 SL combo we'll use to make a trial helix
00101   for (int iSegCombo = 0; iSegCombo < MdcxParameters::nSegCombo; iSegCombo++) {
00102     //FIXME optimize: group segList by superlayer ahead
00103     int axlay = MdcxParameters::findTrkGroup[iSegCombo][0];
00104     int stUlay = MdcxParameters::findTrkGroup[iSegCombo][1];
00105     int stVlay = MdcxParameters::findTrkGroup[iSegCombo][2];
00106 
00107     if(4 == m_debug) std::cout <<" Test combo: ("<< axlay<<","<<stUlay<<","<<stVlay <<")"<< std::endl; 
00108 
00109     float max_dPhiAU = MdcxParameters::maxDp12[iSegCombo];
00110     float max_dPhiAV = MdcxParameters::maxDp13[iSegCombo];
00111 
00112     //---------------------------------------------------------------------
00113     //  Axial SL segment furnishes initial guess for d0, phi0, omega
00114     //---------------------------------------------------------------------
00115     for (int iAx = 0; iAx < nSeg; iAx++) {
00116       bool b_useGood_stU = true; //yzhang add 2009-10-21 17:20:45
00117       bool b_useGood_stV = true; //yzhang add 2009-10-21 17:20:45
00118       if ((segList[iAx]->GetUsedOnHel() != 0) || (segList[iAx]->Fail() != 0) 
00119           || (segList[iAx]->SuperLayer(0) != axlay) || (segList[iAx]->Origin() != -1) ) continue;
00120 
00121       if (4 == m_debug){
00122         std::cout<< "---------1.Ax seg------ No."<<iAx<< std::endl;
00123         segList[iAx]->printSegAll();
00124       }
00125 
00126       //Skip low pt track
00127       omegag = segList[iAx]->Omega();
00128       if (fabs(omegag) > MdcxParameters::maxTrkOmega) {//FIXME 0.2
00129         if (4 == m_debug) std::cout<< "SKIP by maxTrkOmega "<<fabs(omegag) << std::endl;
00130         continue;  
00131       }
00132       if(4 == m_debug) std::cout << " iAx omegag = " << omegag << std::endl;
00133       if(g_omegag) g_omegag->fill(omegag);
00134 
00135       const HepAList<MdcxHit>& hitsAx = segList[iAx]->XHitList();
00136 
00137       //Interpret segment as straight line also
00138       d0g_sl = segList[iAx]->D0_sl_approx();
00139       phi0g_sl = segList[iAx]->Phi0_sl_approx(); 
00140       omegag_sl = 0.0;
00141 
00142       //save c-tors; do it locally
00143       sinPhi0g_sl = -sin(phi0g_sl);
00144       cosPhi0g_sl = cos(phi0g_sl);
00145       xp_sl = segList[iAx]->Xref() + sinPhi0g_sl*d0g_sl;
00146       yp_sl = segList[iAx]->Yref() + cosPhi0g_sl*d0g_sl;
00147       d0g_sl = yp_sl*cosPhi0g_sl + xp_sl*sinPhi0g_sl;
00148       fltL_sl = segList[iAx]->Xref()*cosPhi0g_sl - segList[iAx]->Yref()*sinPhi0g_sl; 
00149       if(4 == m_debug) {
00150         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; 
00151       }
00152 
00153       if (fabs(omegag) > Constants::epsilon) {
00154         MdcxHel ohel = TakeToOrigin(*segList[iAx]);
00155         omegag = ohel.Omega();
00156         phi0g = ohel.Phi0();
00157         d0g = ohel.D0(); 
00158         xc = ohel.Xc();
00159         yc = ohel.Yc();
00160         x0g = ohel.X0();
00161         y0g = ohel.Y0();
00162         phiAx = atan2(y0g-yc, x0g-xc);
00163         fltAx = dlen(phi0g,segList[iAx]->Phi0(), omegag); 
00164         if(4 == m_debug) {
00165           std::cout <<"--Ax :"<< " ohel ";
00166           ohel.print();
00167         }
00168       }
00169 
00170       //---------------------------------------------------------------------
00171       //First stereo SL segement and axial segment furnish one z intersection
00172       //---------------------------------------------------------------------
00173       //for (int iStU = 0; iStU < nSeg; iStU++ ) //yzhang delete 2009-10-21 18:13:50
00174       for (int iStU = -1; true; ) {//yzhang add
00175         if(!b_useGood_stU && (iStU >= (nSeg-1))) break;
00176         if(b_useGood_stU && (iStU >= (nSeg-1))){
00177           iStU = -1;
00178           b_useGood_stU = false;
00179         }
00180         iStU++;
00181         if ((segList[iAx]->GetUsedOnHel() != 0) 
00182             || (segList[iStU]->GetUsedOnHel() != 0) 
00183             || (segList[iStU]->Fail() != 0) 
00184             || (segList[iStU]->SuperLayer(0) != stUlay) 
00185             || (segList[iStU]->Origin() != -1) 
00186             //yzhang add  2009-10-21 17:21:55
00187             || (b_useGood_stU && (segList[iStU]->Nhits()< 4))
00188            ){ continue;}
00189 
00190         if (4 == m_debug){
00191           std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl; 
00192           std::cout<< "---------2.St U seg ------No."<<iStU<< std::endl;//yzhang debug
00193           segList[iStU]->printSeg();
00194           std::cout<<" omega of seg iStU = "<<segList[iStU]->Omega()<<std::endl;
00195         }
00196 
00197         //Test dPhi between Ax slayer and StU slayer 
00198         const HepAList<MdcxHit>& hitsStU = segList[iStU]->XHitList();
00199         double dPhiAU = fabs(hitsAx[0]->phiMid() - hitsStU[0]->phiMid());
00200         if (dPhiAU > Constants::pi) dPhiAU = Constants::twoPi - dPhiAU;
00201         if(4 == m_debug) {
00202           if (dPhiAU > max_dPhiAU){
00203             cout << "SKIP by dPhiAU " <<dPhiAU <<" > "<< max_dPhiAU << endl;
00204           }else{
00205             cout << "KEEP by dPhiAU " <<dPhiAU<< " <= " << max_dPhiAU << endl;
00206           }
00207           std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
00208             <<" StU "<<hitsStU[0]->phiMid() << std::endl; 
00209           std::cout<<  std::endl;
00210         }
00211         if (g_dPhiAU) g_dPhiAU->fill(dPhiAU,segList[iStU]->SuperLayer());
00212 
00213         if (dPhiAU > max_dPhiAU){ continue; }
00214 
00215         //Calculate z by Ax and StU segs
00216         xl1_0 = segList[iStU]->Xline_bbrrf();
00217         yl1_0 = segList[iStU]->Yline_bbrrf();
00218         sx1   = segList[iStU]->Xline_slope();
00219         sy1   = segList[iStU]->Yline_slope();
00220         if (4 == m_debug) std::cout prt(xl1_0) prt(yl1_0) prt(sx1) prt(sy1) prt(omegag) << std::endl;
00221         if (fabs(omegag) >Constants::epsilon) {
00222           zintAU = findz_cyl(iAx, iStU, segList);
00223           xl1 = xl1_0 + sx1*zintAU;
00224           yl1 = yl1_0 + sy1*zintAU;
00225           rl1 = sqrt(xl1*xl1 + yl1*yl1);
00226           phiStU = atan2(yl1-yc, xl1-xc);
00227           if (4 == m_debug) cout prt(zintAU) prt(xl1) prt(yl1) prt(rl1) prt(phiStU) <<endl;
00228         } else {
00229           zintAU = -9999.;
00230         }
00231 
00232         //Calculate z_sl by Ax and StU segs
00233         zintAU_sl = findz_sl(iAx, iStU, segList); 
00234         xl1_sl = xl1_0 + sx1*zintAU_sl;
00235         yl1_sl = yl1_0 + sy1*zintAU_sl;
00236         rl1_sl = sqrt(xl1_sl*xl1_sl + yl1_sl*yl1_sl);
00237         if (4 == m_debug) cout prt(zintAU_sl) prt(xl1_sl) prt(yl1_sl) prt(rl1_sl) <<endl;
00238 
00239         //Cut by z and z_sl
00240         if ( (zintAU < -MdcxParameters::maxMdcZLen) &&
00241             (zintAU_sl < -MdcxParameters::maxMdcZLen) ) {
00242           if (4 == m_debug) std::cout << " SKIP by zintAU < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
00243           continue;
00244         }
00245         if ( (zintAU > MdcxParameters::maxMdcZLen) &&
00246             (zintAU_sl > MdcxParameters::maxMdcZLen) ) {
00247           if (4 == m_debug) std::cout << " SKIP by zintAU > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
00248           continue;
00249         }
00250 
00251         //---------------------------------------------------------------------
00252         //Next stereo SL segement and axial segment furnish other z intersection
00253         //---------------------------------------------------------------------
00254         for (int iStV = -1; true; ) {
00255           if(!b_useGood_stV && (iStV >= (nSeg-1))) break;
00256           if(b_useGood_stV && (iStV >= (nSeg-1))){
00257             iStV = -1;
00258             b_useGood_stV = false;
00259           }
00260           iStV++;
00261           if ((segList[iStV]->GetUsedOnHel()!=0) 
00262               || (segList[iStU]->GetUsedOnHel()!=0) 
00263               || (segList[iAx]->GetUsedOnHel()!=0) 
00264               || (segList[iStV]->Fail() != 0) 
00265               || (segList[iStV]->SuperLayer(0) != stVlay) 
00266               || (segList[iStV]->Origin() != -1)  
00267               //yzhang add  2009-10-21 17:21:55
00268               || (b_useGood_stV && (segList[iStV]->Nhits()< 4))
00269              ){ continue; }
00270 
00271           if (4 == m_debug){
00272             std::cout <<" Test combo: AUV "<< axlay<<","<<stUlay<<","<<stVlay << std::endl; 
00273             std::cout<< "---------3.St V seg ------No."<<iStV<< std::endl;//yzhang debug
00274             segList[iStV]->printSeg();
00275             std::cout<<" omega of seg iStV = "<<segList[iStV]->Omega()<<std::endl;
00276           }
00277 
00278           //Calculate dPhi between Ax and StV segs
00279           const HepAList<MdcxHit>& hitsStV = segList[iStV]->XHitList();
00280           double dPhiAV = fabs(hitsAx[0]->phiMid() - hitsStV[0]->phiMid());
00281           if (dPhiAV > Constants::pi) dPhiAV = Constants::twoPi - dPhiAV;
00282 
00283           if(4 == m_debug) {
00284             if (dPhiAV > max_dPhiAV){
00285               cout << "SKIP by dPhiAV " <<dPhiAV <<" > "<< max_dPhiAV << endl;
00286             }else{
00287               cout << "KEEP by dPhiAV " <<dPhiAV<< " <= " << max_dPhiAV << endl;
00288             }
00289             std::cout << " phi mid Ax "<<hitsAx[0]->phiMid()
00290               <<" StV "<<hitsStV[0]->phiMid() << std::endl; 
00291             std::cout<<  std::endl;
00292           }
00293           if (g_dPhiAV) g_dPhiAV->fill(dPhiAV,segList[iStV]->SuperLayer());
00294           if (dPhiAV > max_dPhiAV) { continue; }
00295 
00296 
00297           //Calculate z by Ax and StV segs
00298           xl2_0 = segList[iStV]->Xline_bbrrf();
00299           yl2_0 = segList[iStV]->Yline_bbrrf();
00300           sx2   = segList[iStV]->Xline_slope();
00301           sy2   = segList[iStV]->Yline_slope();
00302           if (fabs(omegag) > Constants::epsilon) {
00303             zintAV = findz_cyl(iAx, iStV, segList);
00304             xl2 = xl2_0 + sx2*zintAV;
00305             yl2 = yl2_0 + sy2*zintAV;
00306             rl2 = sqrt(xl2*xl2 + yl2*yl2);
00307             if (4 == m_debug){ 
00308               segList[iAx]->printSeg();
00309               segList[iStV]->printSeg();
00310               cout << "zintAV " << zintAV << endl;
00311             }
00312             phiStV = atan2(yl2-yc, xl2-xc);
00313           } else {
00314             zintAV = -9999.;
00315           }
00316 
00317           //Calculate z by Ax and StV segs
00318           zintAV_sl = findz_sl(iAx, iStV, segList);
00319           xl2_sl = xl2_0 + sx2*zintAV_sl;
00320           yl2_sl = yl2_0 + sy2*zintAV_sl; 
00321           rl2_sl = sqrt(xl2_sl*xl2_sl + yl2_sl*yl2_sl);
00322 
00323           //Cut by z of AV segs
00324           if ( (zintAV < -MdcxParameters::maxMdcZLen) &&
00325               (zintAV_sl < -MdcxParameters::maxMdcZLen) ) {
00326             if (4 == m_debug) std::cout << " SKIP by zintAV < max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
00327             continue;
00328           }
00329           if ( (zintAV > MdcxParameters::maxMdcZLen) &&
00330               (zintAV_sl > MdcxParameters::maxMdcZLen) ) {
00331             if (4 == m_debug) std::cout << " SKIP by zintAV > max MdcZLen "<<MdcxParameters::maxMdcZLen << std::endl; 
00332             continue;
00333           }
00334 
00335           //We now have enough info to form and fit a guess helix (inside iAx-iStU-iStV loops)
00336           if (4 == m_debug) {
00337             //cout << "KEEP layer[klm]: [" << axlay <<","<< stUlay ","<< stVlay << "] " prt(d0g) prt(phi0g) prt(omegag) zrt(zintAU) prt(zintAV) << endl; 
00338             segList[iAx]->printSeg();
00339             segList[iStU]->printSeg();
00340             segList[iStV]->printSeg();
00341           }
00342           MdcxFittedHel fithel;
00343           double helixFitSigma = MdcxParameters::helixFitSigma;
00344           double first_prob = 0.0; 
00345           HepAList<MdcxHit> hitlist;
00346           hitlist.append(hitsAx);
00347           hitlist.append(hitsStU);
00348           hitlist.append(hitsStV);
00349           if(g_trkllmk) for (int i=0; i<hitlist.length(); i++ ){ g_trkllmk->fill(hitlist[i]->Layer()); }
00350 
00351 
00352           //---------------------------------------------------------------------
00353           //Try fitting the trial helix from axial SL as circle first
00354           //---------------------------------------------------------------------
00355           floop = 1;
00356           while (floop) {
00357             if(4 == m_debug) std::cout<< "---------4.Ax circle fit------"<< std::endl;//yzhang debug
00358             if (fabs(omegag) < Constants::epsilon) {
00359               if(4 == m_debug) std::cout<<"SKIP by omegag==0   "<<std::endl;
00360               break;
00361             }
00362             if ( (zintAU < -MdcxParameters::maxMdcZLen) ||
00363                 (zintAU > MdcxParameters::maxMdcZLen) ){
00364               if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
00365               break; 
00366             }
00367             if ( (zintAV < -MdcxParameters::maxMdcZLen) ||
00368                 (zintAV > MdcxParameters::maxMdcZLen) ){
00369               if(4 == m_debug) std::cout<<"SKIP by zintAU out of range "<<zintAU<<" "<<MdcxParameters::maxMdcZLen<<std::endl;
00370               break; 
00371             }
00372 
00373 
00374             //Calculate flight length of AUV segs
00375             double fltLenUV = dlen(phiStU, phiStV, omegag);
00376             if(4 == m_debug){
00377               std::cout << "MdcxFindTracks "  prt(fltLenUV) prt(phiStU) prt(phiStV) prt(omegag)<< std::endl;
00378             }
00379             if (fltLenUV > MdcxParameters::maxDlen) {//15. to 150.
00380               if(4 == m_debug) std::cout<<"SKIP by fltLenUV > "<<MdcxParameters::maxDlen<<std::endl;
00381               break;   //FIXME
00382             }
00383             tanlg = (zintAV - zintAU)/fltLenUV; 
00384             double fltLenAU = dlen(phiAx, phiStU, omegag);
00385             z0g = zintAU - fltLenAU*tanlg;
00386             if (4 == m_debug) cout prt(z0g) prt(tanlg) prt(fltLenUV) prt(zintAU) prt(zintAV)<< endl;
00387 
00388             //Cut by z0g
00389             if ((z0g < -MdcxParameters::maxMdcZLen) || (z0g > MdcxParameters::maxMdcZLen)) {
00390               if (4 == m_debug) std::cout<<"SKIP by z0g out of range "<<z0g<<">"<<MdcxParameters::maxMdcZLen<<std::endl;
00391               break;
00392             }
00393             ztest = z0g + fltAx*tanlg; 
00394             if (4 == m_debug) std::cout prt(ztest) prt(fltAx)<<std::endl;
00395 
00396             //---Helix fitting
00397             //MdcxHel ghel(d0g,phi0g,omegag,z0g,tanlg); // ghel.print();
00398             MdcxHel ghel(segList[iAx]->D0(), segList[iAx]->Phi0(), segList[iAx]->Omega(),
00399                 ztest, tanlg, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref()); 
00400             //zoujh?: ghel.SetTurnFlag(1); // ghel.print();
00401             MdcxFittedHel firstfit(hitlist,ghel,helixFitSigma);  
00402             first_prob = firstfit.Prob(); 
00403             if (4 == m_debug) {
00404               std::cout<<" after first fit:  "<<std::endl;
00405               firstfit.FitPrint();
00406             }
00407             if (first_prob >= probmin) {
00408               MdcxHel helixOrigin= TakeToOrigin(firstfit);
00409               MdcxFittedHel nextfit(hitlist,helixOrigin,helixFitSigma);
00410               first_prob = nextfit.Prob();
00411               fithel=nextfit;
00412               if (g_trklgood) for (int i=0; i<nextfit.Nhits(); i++){ g_trklgood->fill(nextfit.Layer(i));}
00413             }
00414             if ( g_trklfirstProb) g_trklfirstProb->fill(first_prob);
00415             if (4 == m_debug) cout << " prob after helix fitting = " << first_prob << endl;     
00416             floop = 0;
00417           }//end fitting trial helix
00418 
00419           //---------------------------------------------------------------------
00420           // Try fitting straight line trial helix if fit from axial-SL-as-circle trial
00421           //  if helix not good enough
00422           //---------------------------------------------------------------------
00423           floop = 1;
00424           while (floop) {
00425             if (4 == m_debug) std::cout<< "---------4.2 straight line fit------"<< std::endl;//yzhang debug
00426             if (4 == m_debug) cout << " zintAU_sl " << zintAU_sl << " zintAV_sl " << zintAV_sl << endl; 
00427             if (first_prob > probmin) {
00428               if (4 == m_debug) cout << " helix fit good , exit straight line fit." << endl;     
00429               break;
00430             }
00431 
00432             if ( (zintAU_sl < -MdcxParameters::maxMdcZLen) ||
00433                 (zintAU_sl > MdcxParameters::maxMdcZLen) ) break;
00434             if ( (zintAV_sl < -MdcxParameters::maxMdcZLen) ||
00435                 (zintAV_sl > MdcxParameters::maxMdcZLen) ) break;
00436 
00437             double dx = xl2_sl - xl1_sl;
00438             double dy = yl2_sl - yl1_sl;
00439             double dl = sqrt(dx*dx + dy*dy);
00440             tanlg_sl = (zintAV_sl - zintAU_sl)/dl;
00441             dx = xl1_sl + d0g_sl*sin(phi0g_sl);
00442             dy = yl1_sl - d0g_sl*cos(phi0g_sl);
00443             dl = sqrt(dx*dx + dy*dy);
00444             z0g_sl = zintAU_sl - dl*tanlg_sl;
00445             if (4 == m_debug) cout prt(d0g_sl) prt(phi0g_sl) prt(z0g_sl) prt(tanlg_sl)<< endl;
00446 
00447             if ((z0g_sl < -MdcxParameters::maxMdcZLen) || (z0g_sl > MdcxParameters::maxMdcZLen)){
00448               if(4 == m_debug) std::cout << "SKIP by z0g_sl:"<<z0g_sl<<" > "<<MdcxParameters::maxMdcZLen << std::endl; 
00449               break;
00450             }
00451             //MdcxHel ghel(d0g_sl,phi0g_sl,omegag_sl,z0g_sl,tanlg_sl); // ghel.print();
00452             ztest_sl = z0g_sl + fltL_sl*tanlg_sl;
00453             if (4 == m_debug) cout prt(ztest_sl) << endl;
00454 
00455             //Helix fitting
00456             MdcxHel ghel(segList[iAx]->D0_sl_approx(), segList[iAx]->Phi0_sl_approx(), omegag_sl,
00457                 ztest_sl, tanlg_sl, 0.0, 11111, 0, segList[iAx]->Xref(), segList[iAx]->Yref());
00458             //zoujh?: ghel.SetTurnFlag(1); // ghel.print();
00459             MdcxFittedHel firstfit(hitlist, ghel, helixFitSigma);
00460             first_prob = firstfit.Prob();
00461             //if (first_prob >= probmin)fithel=firstfit;
00462             if(4 == m_debug){
00463               ghel.print();
00464               std::cout <<"first_prob "<<first_prob << std::endl;  
00465             }
00466 
00467             if (first_prob >= probmin) {
00468               MdcxHel helixOrigin = TakeToOrigin(firstfit);
00469               MdcxFittedHel nextfit(hitlist, helixOrigin, helixFitSigma);
00470               first_prob = nextfit.Prob();
00471               fithel = nextfit;
00472               if (4 == m_debug) {
00473                 cout << " three seg sl nextfit" <<  endl;
00474                 nextfit.FitPrint();
00475               }
00476             }
00477             floop = 0;
00478           }
00479 
00480           //-----------------------------------------------------------------
00481           //Got a good trial helix fit; now try to add other segments onto it
00482           //-----------------------------------------------------------------
00483           floop = 1;
00484           while (floop) {
00485             floop = 0;
00486             if(4 == m_debug)std::cout<< "---------5. try add seg to helix------"<< std::endl;
00487             if (first_prob < probmin) {
00488               if(4 == m_debug) std::cout<< "prob"<<first_prob<< " "
00489                 <<probmin<<" fit NOT good, exit add segs "<< std::endl;
00490               break;
00491             }
00492             float bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
00493             int   bndof = fithel.Nhits() - 5;
00494             float bprob = Mdcxprobab(bndof, bchisq);
00495             trkCount++; 
00496             segList[iAx]->SetUsedOnHel(trkCount);
00497             segList[iStU]->SetUsedOnHel(trkCount);
00498             segList[iStV]->SetUsedOnHel(trkCount); 
00499 
00500             if (4 == m_debug){
00501               cout << "in add seg to trail helix, klm seg:" << endl;  
00502               segList[iAx]->printSeg();                  
00503               segList[iStU]->printSeg();                  
00504               segList[iStV]->printSeg();                  
00505             }
00506 
00507             //Loop superlayers
00508             for (int iSlay = 0; iSlay < Constants::nSuperLayer; iSlay++) {
00509               int ilay = MdcxParameters::layerSet2AddSeg[iSegCombo][iSlay]; 
00510               for (int iSeg = 0; iSeg < nSeg; iSeg++) {
00511                 if ((segList[iSeg]->SuperLayer(0) != ilay) 
00512                     ||(segList[iSeg]->GetUsedOnHel() != 0)
00513                     || (segList[iSeg]->Fail() != 0) 
00514                     || (segList[iSeg]->Origin() != -1)) continue;
00515                 if(4 == m_debug) {
00516                   std::cout<< "try to add seg No."<< iSeg<<" to track  "<< std::endl;
00517                   segList[iSeg]->printSeg();
00518                 }
00519 
00520                 //Cut by phi_diff = (Phi_Ax - Phi_Add)
00521                 float phiAxSeg = segList[iAx]->XHitList()[0]->phiMid();
00522                 float phiAddSeg = segList[iSeg]->XHitList()[0]->phiMid();
00523                 double phiKNSegDiff = fabs(phiAxSeg - phiAddSeg);
00524                 if (phiKNSegDiff > 0.5*Constants::pi && phiKNSegDiff < 1.5*Constants::pi ) {
00525                   if(4 == m_debug)std::cout<< " SKIP by seg phi diff of phi ,iAx:"<<
00526                     phiAxSeg<<" iSeg:"<<phiAddSeg<<" diff:"<<phiKNSegDiff << std::endl;//yzhang debug 
00527                   continue;  // zoujh FIXME
00528                 }
00529 
00530                 //Calculate doca
00531                 xp = segList[iSeg]->Xref() - segList[iSeg]->SinPhi0()*segList[iSeg]->D0();
00532                 yp = segList[iSeg]->Yref() + segList[iSeg]->CosPhi0()*segList[iSeg]->D0();
00533                 const HepAList<MdcxHit>& hitsSegAdd= segList[iSeg]->XHitList();
00534                 double proca = fithel.Doca( hitsSegAdd[0]->wx(), hitsSegAdd[0]->wy(),
00535                     hitsSegAdd[0]->wz(), xp, yp );
00536                 if (g_trklproca) g_trklproca->fill(proca);
00537 
00538                 //Segment passes loose doca cut; does it belong on this track?
00539                 if (!((fabs(proca) < MdcxParameters::maxProca) &&(fithel.Doca_Len() < fithel.Lmax())) ){
00540                   //if (g_trklprocaSl) g_trklprocaSl->fill(segList[iSeg]->SuperLayer(0));
00541                   if(4 == m_debug){
00542                     std::cout<< "SKIP by fabs(proca):"<< fabs(proca)<< " maxProca:"<<MdcxParameters::maxProca <<" Doca_Len:"<<fithel.Doca_Len() <<"<Lmax:"<<fithel.Lmax()<< std::endl;//yzhang debug 
00543                   }
00544                 }
00545 
00546                 if ((fabs(proca) < MdcxParameters::maxProca) && (fithel.Doca_Len() < fithel.Lmax())) {
00547                   MdcxFittedHel newhel;
00548                   newhel.Grow(fithel, hitsSegAdd);
00549                   if (newhel.Nhits() != hitlist.length()) {
00550                     if (4 == m_debug){
00551                       cout << "try to add seg " << iSeg  
00552                         << ": prob " << newhel.Prob()
00553                         <<" rcs "<<newhel.Rcs()<<" <=? "
00554                         <<MdcxParameters::maxRcsInAddSeg << endl;
00555                     }
00556                     //yzhang change 2009-10-21 FIXME
00557                     //if (newhel.Rcs() > MdcxParameters::maxRcsInAddSeg)
00558                     if (newhel.Prob() < probmin){
00559                       if(4 == m_debug){ cout << "helix.prob " << newhel.Prob() 
00560                         << "<? "<<probmin<<"  Fail=" << newhel.Fail() 
00561                           <<". and Call ReFit()"<< endl;  
00562                       }
00563                       newhel.ReFit(); //FIXME
00564                     }
00565                     if(4 == m_debug) {
00566                       cout << "newhel.prob " << newhel.Prob() 
00567                         << "<probmin "<<probmin<<" Fail=" << newhel.Fail() << endl;
00568                     }
00569                     //yzhang change 2009-10-21 FIXME
00570                     if ( (newhel.Prob() >= probmin) && (newhel.Fail() == 0) ) {
00571                     //if ( (newhel.Rcs() <= MdcxParameters::maxRcsInAddSeg) && (newhel.Fail() == 0) ) 
00572                       fithel = newhel;
00573                       hitlist = newhel.XHitList();
00574                       segList[iSeg]->SetUsedOnHel(trkCount);
00575                       bchisq = fithel.Chisq()*helixFitSigma*helixFitSigma;
00576                       bndof = fithel.Nhits() - 5;
00577                       float tprob = Mdcxprobab(bndof, bchisq); 
00578                       if (tprob > bprob) bprob = tprob;
00579                       if (4 == m_debug) {
00580                         cout << "add seg " << iSeg << " with new hits; prob now " << newhel.Prob() << endl;
00581                         cout << "  add "; segList[iSeg]->printSeg();  
00582                       }
00583                     }   //*end trial track still good so keep new seg 
00584                     else{
00585                       if(4 == m_debug)std::cout<<" helix fit failed"<<std::endl;
00586                     }
00587                   } else { //*end new hits so refit
00588                     segList[iSeg]->SetUsedOnHel(trkCount);
00589                     if (4 == m_debug) cout << "segment " << iSeg << " added but no new hits" << endl;
00590                   }//*end dont need to refit if no new hits
00591                 }//*end fabs(proca) < 0.010
00592               }//*end iSeg (loop over segs in a superlayer)
00593             }//*end iSlay (loop over superlayers)
00594             if (((fithel.Prob() < 0.90) && (fithel.Nhits() == 12)) || (fithel.Fail() != 0)) {
00595               for (int i = 0; i < nSeg; i++) { if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0); }
00596               trkCount--;
00597               break; //*jump to end floop first_prob >= probmin
00598             }
00599 
00600             if(4 == m_debug){
00601               std::cout<< "After add segs" << std::endl;//yzhang debug
00602               //fithel.printAll();
00603             }
00604 
00605             //See if helix better moving backwards
00606             if(4 == m_debug) std::cout<< "---------6. flip------"<< std::endl;//yzhang debug
00607             fithel.flip();
00608 
00609             //See if this track shares lots of hits with another saved track
00610             int kki = 0;
00611             int notduplicate = 1; 
00612             while(MdcxTrklist[kki]) {
00613               if (4 == m_debug)std::cout<< "---------7. test saved track------"<< std::endl;//yzhang debug
00614               const HepAList<MdcxHit>& junk=MdcxTrklist[kki]->XHitList();
00615               int overl = 0; 
00616               int kkj = 0;
00617               while (junk[kkj]) {
00618                 int kkl = 0;
00619                 while (hitlist[kkl]) {
00620                   if (hitlist[kkl] == junk[kkj]) overl++; 
00621                   kkl++;
00622                 }
00623                 kkj++;
00624               }
00625               if (4 == m_debug) cout << "overlap with track # " << kki << " is "
00626                 << junk.length() << " hitList " << hitlist.length() << " overlap " << overl << endl;
00627               if ( (hitlist.length()/4.) <= overl) notduplicate = 0;
00628               //yzhang 2009-10-21 16:30:53
00629               //if ( hitlist.length() <= (overl+2) ) notduplicate = 0;
00630               kki++;
00631             }
00632 
00633             if(g_trklhelix) for (int iHit=0; iHit<fithel.Nhits(); iHit++){ g_trklhelix->fill(fithel.Layer(iHit)); }
00634 
00635             //If not a duplicate, is it worth saving?
00636             if (notduplicate) {
00637               if(4 == m_debug) std::cout<< "---------8. test worth saving?------"<< std::endl;//yzhang debug
00638               MdcxFittedHel *finehel = new MdcxFittedHel(hitlist, fithel);
00639 
00640               if(4 == m_debug) {
00641                 std::cout<<__FILE__
00642                 <<" finehel Prob "<<finehel->Prob()
00643                   <<" nhits "<<finehel->Nhits()
00644                   <<" bprob "<< bprob 
00645                   <<" probmin "<< probmin
00646                   <<" Fail "<< finehel->Fail()
00647                   <<  std::endl;
00648                 finehel->FitPrint();
00649               }
00650               //cut nhits from 15->9 yzhang 2009-10-21 
00651               if (( (finehel->Prob()>0.001) || (finehel->Nhits()>=15) || (bprob>probmin) 
00652                     //|| (bchisq/bndof > MdcxParameters::maxRcsInAddSeg) 
00653                   ) && (finehel->Fail() == 0) ) {
00654                 nfound++;
00655                 sumprob += finehel->Prob();
00656                 //if ( hitlist.length() == 16 )resout(finehel);
00657 
00658                 //Begin hit dropping section
00659                 //Set nodrop != 0 if want to prohibit dropping bad hits 
00660                 int nodrop = 0; 
00661                 if ((finehel->Prob() > 0.05) || nodrop ) {
00662                   MdcxTrklist.append(finehel);
00663                   //yzhang hist
00664                   for (int iHit=0; iHit<finehel->Nhits(); iHit++){
00665                     if(g_trklappend1) g_trklappend1->fill(finehel->Layer(iHit));//yzhang hist
00666                   }
00667                   //zhangy hist
00668                 } else {
00669                   MdcxFittedHel* drophel = new MdcxFittedHel( drophits(finehel) );
00670                   float ratdrop = float(drophel->Nhits()) / float(finehel->Nhits());
00671                   if ((drophel->Prob() > finehel->Prob()) && 
00672                       (ratdrop > 0.50) && (drophel->Fail() == 0)) {  //FIXME
00673                     if (4 == m_debug) cout << "Track " << trkCount << " drops " 
00674                       << finehel->Nhits()-drophel->Nhits() << " hits" << endl;
00675                     MdcxTrklist.append(drophel);
00676                     if(g_trklappend2) for (int iHit=0;iHit<drophel->Nhits();iHit++){g_trklappend2->fill(drophel->Layer(iHit));}
00677                     delete finehel;
00678                   } else {
00679                     MdcxTrklist.append(finehel);
00680                     if(g_trklappend3)for(int iHit=0;iHit<finehel->Nhits();iHit++){g_trklappend3->fill(finehel->Layer(iHit));}
00681                     delete drophel;
00682                   }
00683                 }
00684 
00685                 //End hit dropping section
00686                 int nwl = MdcxTrklist.length() - 1;
00687                 if ((4 == m_debug) && (nwl > -1)) {
00688                   if (4 == m_debug) {
00689                     cout << endl << "found track " << trkCount<<std::endl;
00690                     MdcxTrklist[nwl]->print();
00691                     MdcxTrklist[nwl]->FitPrint();
00692                   }
00693                 }
00694               } else { //*endif good fit so add to MdcxTrklist
00695                 for (int i = 0; i < nSeg; i++) {
00696                   if (segList[i]->GetUsedOnHel() == trkCount) segList[i]->SetUsedOnHel(0);
00697                 }
00698                 delete finehel;
00699                 trkCount--;
00700               }
00701             }//*end if not duplicate 
00702           }//*end floop first_prob >= probmin
00703         }// end iStV
00704       }// end iStU  
00705     }// end iAx
00706     if(4==m_debug) cout<< "end of this combo"<<endl;
00707   }// end iSegCombo
00708   if(4== m_debug) cout << " In MdcxFindTracks, found " << trkCount << " good tracks" << endl;
00709 }//end of process

void MdcxFindTracks::resout (  MdcxFittedHel *  f  )  [protected]
 

void MdcxFindTracks::resout (  MdcxFittedHel *  f  )  [protected]
  

00828                                                   {
00829   MdcxHel mdcxHel = (MdcxHel) *finehel;
00830   HepAList<MdcxHit>& hitlist = (HepAList<MdcxHit>&)finehel->XHitList();
00831   int kkk = 0; while(hitlist[kkk]){hitlist[kkk]->SetConstErr(0); kkk++;}
00832   MdcxFittedHel newhel(hitlist, mdcxHel);
00833   if(4 == m_debug) newhel.FitPrint(mhel,cout);
00834   kkk = 0;
00835   while(hitlist[kkk]) {
00836     float doca = newhel.Doca(*hitlist[kkk]);
00837     float zh   = newhel.Doca_Zh();
00838     float tof  = newhel.Doca_Tof();
00839     float dd   = hitlist[kkk]->d(newhel);
00840     float tt   = hitlist[kkk]->tcor(zh, tof);  
00841     if (4 == m_debug){
00842       cout << "MdcxFindTracks::resout (" 
00843         << hitlist[kkk]->Layer() << "," 
00844         << hitlist[kkk]->WireNo() << ") dt " 
00845         << tt << " resid " 
00846         << fabs(dd)-fabs(doca) << " doca " 
00847         << fabs(doca) << endl;
00848     }
00849     kkk++;
00850   }
00851   kkk=0; 
00852   while(hitlist[kkk]){hitlist[kkk]->SetConstErr(1); kkk++;}
00853 }//endof drophits

MdcxHel MdcxFindTracks::TakeToOrigin (  MdcxHel &   )  [protected]
 

MdcxHel MdcxFindTracks::TakeToOrigin (  MdcxHel &   )  [protected]
  

00855                                                   {
00856   double lng;
00857   double omegag = ihel.Omega();
00858   double phi0g = ihel.Phi0();
00859   double phi0s = phi0g; 
00860   double xp = ihel.Xref() - ihel.SinPhi0()*ihel.D0();
00861   double yp = ihel.Yref() + ihel.CosPhi0()*ihel.D0();
00862   double d0g = yp*cos(phi0g) - xp*sin(phi0g);
00863   if (omegag != 0.0) {
00864     phi0g = phi0g - Constants::pi;
00865     double xc = sin(phi0g)/omegag;
00866     double yc = -cos(phi0g)/omegag;
00867     double t1 = -xc - xp;
00868     double t2 = yc + yp;
00869     double phinew = atan2(t1, t2);
00870     if (omegag > 0.0) phinew += Constants::pi;
00871     if (phinew > Constants::pi) phinew -= Constants::twoPi;
00872     double xh = xc - sin(phinew)/omegag;
00873     double yh = yc + cos(phinew)/omegag;
00874     double x0 = xh + xp;
00875     double y0 = yh + yp;
00876     d0g = sqrt(x0*x0 + y0*y0);
00877     phi0g = phinew + Constants::pi;  
00878     double sd0 = x0*sin(phi0g) - y0*cos(phi0g);
00879     if (sd0 > 0.0) d0g = -d0g;
00880     lng = dlen(phi0g, phi0s, omegag); 
00881   } else {
00882     lng = ihel.Xref()*ihel.CosPhi0() + ihel.Yref()*ihel.SinPhi0(); 
00883   }
00884   double tanlg = ihel.Tanl();
00885   double z0g = ihel.Z0() - lng*tanlg; 
00886   //cout << " z0g, tanlg, lng " << z0g << " " << tanlg << " " << lng << endl;
00887   double t0g = ihel.T0();
00888   int codeg = ihel.Code(); 
00889   MdcxHel ohel(d0g, phi0g, omegag, z0g, tanlg, t0g, codeg);
00890   return ohel;
00891 }

Member Data Documentation

float MdcxFindTracks::curvmax [protected]
 

int MdcxFindTracks::m_debug [protected]
 

bool MdcxFindTracks::m_doSag [protected]
 

HepAList<MdcxFittedHel> MdcxFindTracks::MdcxTrklist [protected]
 

HepAList<MdcxFittedHel> MdcxFindTracks::MdcxTrklist [protected]
 

MdcxHel MdcxFindTracks::mhel [protected]
 

int MdcxFindTracks::nbad [protected]
 

int MdcxFindTracks::nfail [protected]
 

int MdcxFindTracks::ngood [protected]
 

int MdcxFindTracks::nhitsgd [protected]
 

int MdcxFindTracks::nhitsmc [protected]
 

int MdcxFindTracks::npure [protected]
 

float MdcxFindTracks::probmin [protected]
 

The documentation for this class was generated from the following files:
Generated on Wed Feb 2 16:30:13 2011 for BOSS6.5.5 by  doxygen 1.3.9.1