#include <MdcxAddHits.h>
Public Member Functions | |
MdcxAddHits (HepAList< MdcxFittedHel > &f, const HepAList< MdcxHit > &h, float addit=1.5) | |
MdcxAddHits (HepAList< MdcxHel > &f, const HepAList< MdcxHit > &h, float addit=1.5) | |
virtual | ~MdcxAddHits () |
const HepAList< MdcxHit > & | GetAssociates (int i=0) |
Protected Attributes | |
float | addcut |
int | ncalc |
int | nadded |
double | sumpull |
int | kkk |
int | kkl |
int | thot |
int | tuot |
int | ntracks |
int | nhits |
HepAList< MdcxHel > | tttt |
HepAList< MdcxHit > | hhhh |
HepAList< MdcxHit > | listoasses |
Definition at line 33 of file MdcxAddHits.h.
MdcxAddHits::MdcxAddHits | ( | HepAList< MdcxFittedHel > & | f, | |
const HepAList< MdcxHit > & | h, | |||
float | addit = 1.5 | |||
) |
Definition at line 40 of file MdcxAddHits.cxx.
References addcut, hhhh, kkk, kkl, nhits, ntracks, thot, and tttt.
00041 : 00042 ncalc(0), nadded(0), sumpull(0.0), thot(0), tuot(0) 00043 { 00044 ntracks = trkl.length(); 00045 nhits = hitl.length(); 00046 addcut = addit; 00047 // cout << "MdcxAddHits called with " << ntracks << " tracks, " 00048 // << nhits << " hits, addit = " << addit << endl; 00049 if ( (ntracks < 1) || (nhits < 1) ) return; 00050 kkk = 0; 00051 kkl = 0; 00052 while(hitl[kkl]) { 00053 hitl[kkl]->SetUsedOnHel(0); 00054 hhhh.append(hitl[kkl++]); 00055 } 00056 while(trkl[kkk]) { 00057 MdcxHel* thel = trkl[kkk]; 00058 tttt.append(thel); 00059 const HepAList<MdcxHit>& dchitlist = trkl[kkk++]->XHitList(); 00060 thot += dchitlist.length(); 00061 kkl = 0; 00062 while (dchitlist[kkl]) dchitlist[kkl++]->SetUsedOnHel(1); 00063 } 00064 }
MdcxAddHits::MdcxAddHits | ( | HepAList< MdcxHel > & | f, | |
const HepAList< MdcxHit > & | h, | |||
float | addit = 1.5 | |||
) |
Definition at line 66 of file MdcxAddHits.cxx.
References addcut, hhhh, nhits, ntracks, and tttt.
00066 : 00067 ncalc(0), nadded(0), sumpull(0.0) 00068 { 00069 // 00070 // c-tor designed to work with MdcxHitAdder; all MdcxHit's in hitl are 00071 // assumed to be fresh and fair game. 00072 // 00073 ntracks = trkl.length(); 00074 nhits = hitl.length(); 00075 addcut = addit; 00076 // cout << "MdcxAddHits called with " << ntracks << " tracks, " 00077 // << nhits << " hits, addit = " << addit << endl; 00078 if ( (ntracks < 1) || (nhits < 1) ) return; 00079 hhhh = hitl; 00080 tttt = trkl; 00081 }
MdcxAddHits::~MdcxAddHits | ( | ) | [virtual] |
const HepAList< MdcxHit > & MdcxAddHits::GetAssociates | ( | int | i = 0 |
) |
Definition at line 86 of file MdcxAddHits.cxx.
References MdcxHel::D0(), MdcxParameters::debug, MdcxHel::Doca_Len(), MdcxParameters::firstMdcAxialRadius, g_addHitCut, g_addHitCut2d, MdcxHit::GetUsedOnHel(), hhhh, kkl, MdcxHit::Layer(), listoasses, MdcxHel::Lmax(), M_PI, MdcxParameters::maxMdcRadius, MdcxParameters::maxTrkLength, nadded, ncalc, MdcxParameters::nSigAddHitTrk, ntracks, MdcxHel::Omega(), MdcxHel::Ominfl(), MdcxHel::Phi0(), MdcxHit::print(), MdcxHel::print(), MdcxHit::pull(), MdcxHit::pw(), MdcxHit::SetConstErr(), sumpull, tttt, MdcxHel::Xc(), MdcxHel::Xh(), MdcxHel::Yc(), and MdcxHel::Yh().
Referenced by MdcxTrackFinder::FitMdcxTrack().
00086 { 00087 int debug = MdcxParameters::debug; 00088 00089 listoasses.removeAll(); 00090 if ((whichtrack >= ntracks) || (whichtrack < 0)) { 00091 cout << "asked for associates of track " << whichtrack 00092 << " while there are only " << ntracks << " tracks in list " << endl; 00093 return listoasses; 00094 } 00095 00096 double m_2pi = 2.0*M_PI; 00097 MdcxHel* temphel = tttt[whichtrack]; 00098 00099 if(5 == debug) temphel->print(); 00100 00101 double lmax = temphel->Lmax(); 00102 00103 // calc. phim, phip 00104 double gvin = MdcxParameters::firstMdcAxialRadius;//yzhang 2010-05-05 00105 double gvout = MdcxParameters::maxMdcRadius; 00106 double phiex = 0.38; 00107 double phim = -M_PI; 00108 double phip = M_PI; 00109 double rmin = fabs(temphel->D0()); 00110 double rmax, pc, rc, rhel; 00111 if (temphel->Ominfl()) { 00112 rmax = fabs(temphel->D0() + 2.0/temphel->Omega()); 00113 double xc = temphel->Xc(); 00114 double yc = temphel->Yc(); 00115 pc = atan2(yc,xc); 00116 rc = fabs(temphel->D0() + 1.0/temphel->Omega()); 00117 rhel = fabs(1.0/temphel->Omega()); 00118 } else { 00119 rmax = 1000.; 00120 pc = temphel->Phi0(); 00121 rc = 1000.; 00122 rhel = 1000.; 00123 } 00124 if(5 == debug) { 00125 std::cout<< "==Test add hit phi: rmin >?"<<rmin << " gvin "<<gvin 00126 << " rmax >?"<<rmax << " gvout "<<gvout << std::endl; 00127 } 00128 00129 if ((rmin<gvin) && (rmax>gvout)) { 00130 // case A (exiter) 00131 if(5 == debug) std::cout<<" case A (exiter) "<<std::endl; 00132 00133 double len = 0; 00134 double lstep = m_2pi*rhel/16.; 00135 if (lstep>10.) lstep = 10.; 00136 double xl = temphel->Xh(len); 00137 double yl = temphel->Yh(len); 00138 double phil = atan2(yl, xl); 00139 double rl = sqrt(xl*xl + yl*yl); 00140 int nstep = 0; 00141 double philast = phil; 00142 int isin = 0; 00143 int notout = 1; 00144 while ((notout) && (nstep<20)) { 00145 len += lstep; 00146 nstep++; 00147 xl = temphel->Xh(len); 00148 yl = temphel->Yh(len); 00149 phil = atan2(yl, xl); 00150 rl = sqrt(xl*xl + yl*yl); 00151 if ((rl<gvin) && (!isin)) { 00152 philast = phil; 00153 } else if ((rl>gvin) && (!isin)) { 00154 isin = 1; phim = philast; phip = philast; 00155 } 00156 if (isin) { 00157 double deltap = phil - philast; 00158 if (deltap > M_PI) phil -= m_2pi; 00159 if (deltap < -M_PI) phil += m_2pi; 00160 philast = phil; 00161 if (rl > gvout) notout = 0; 00162 if (phil < phim) phim = phil; 00163 if (phil > phip) phip = phil; 00164 } 00165 }//end while 00166 } else if ((rmin>gvin) && (rmax>gvout)) { 00167 if(5 == debug) std::cout<<" case B (albedo) "<<std::endl; 00168 // case B (albedo) 00169 double len=0; double lstep=m_2pi*rhel/16.; if (lstep>10.)lstep=10.; 00170 double xl=temphel->Xh(len); double yl=temphel->Yh(len); 00171 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl); 00172 phim=phil; phip=phil; double phis=phil; double deltap,dp1,dp2; 00173 int nstep=0; double philast=phil; 00174 while ((rl<gvout)&&(nstep<20)){ 00175 len+=lstep; nstep++; xl=temphel->Xh(len); yl=temphel->Yh(len); 00176 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl); 00177 deltap=phil-philast; if (deltap>M_PI)phil-=m_2pi; 00178 if (deltap<-M_PI)phil+=m_2pi; philast=phil; 00179 if (phil<phim)phim=phil; if (phil>phip)phip=phil; 00180 } 00181 dp1=fabs(phis-phim); dp2=fabs(phip-phis); deltap=dp1; 00182 if(dp2>dp1)deltap=dp2; phim=phis-deltap; phip=phis+deltap; 00183 } else if ((rmin<gvin) && (rmax<gvout)) { 00184 // case C (looper) 00185 if(5 == debug) std::cout<<" case C (looper) "<<std::endl; 00186 if (rc>rhel){ 00187 double alp=asin(rhel/rc); phim=pc-alp; phip=pc+alp; 00188 }else{ 00189 // going to have to step to find min, max 00190 double len=0; 00191 double lstep=m_2pi*rhel/16.; 00192 if (lstep>10.)lstep=10.; 00193 if(5 == debug) std::cout<< "ini step "<<m_2pi*rhel/16.<<" lstep " << lstep <<"cm"<<std::endl; 00194 double xl=temphel->Xh(len); double yl=temphel->Yh(len); 00195 double phil=atan2(yl,xl); double rl=sqrt(xl*xl+yl*yl); 00196 int nstep=0; double philast=phil; int isin=0; int notout=1; 00197 while ((notout)&&(nstep<33)){ 00198 len+=lstep; nstep++; xl=temphel->Xh(len); yl=temphel->Yh(len); 00199 phil=atan2(yl,xl); rl=sqrt(xl*xl+yl*yl); 00200 if(5 == debug) { 00201 std::cout<< nstep<<" rl "<<rl<< " gvin " <<gvin<< " isin "<<isin; 00202 } 00203 if ((rl<gvin)&&(!isin)){ 00204 philast=phil; 00205 }else if((rl>gvin)&&(!isin)){ 00206 isin=1; phim=philast; phip=philast; 00207 } 00208 if (isin){ 00209 double deltap=phil-philast; if (deltap>M_PI)phil-=m_2pi; 00210 if (deltap<-M_PI)phil+=m_2pi; philast=phil; if (rl<gvin)notout=0; 00211 if (phil<phim)phim=phil; if (phil>phip)phip=phil; 00212 } 00213 //yzhang add 2011-10-10 00214 if(len > M_PI*rhel*0.75) { 00215 if(5 == debug) { 00216 std::cout<< " len "<<len<<" >pi*R_helix "<<M_PI*rhel<< " rhel "<<rhel<< " break"<<std::endl; 00217 } 00218 break; 00219 } 00220 if(5 == debug) { 00221 std::cout<< " philast "<<philast<<" phim "<<phim << " phip "<<phip <<" len "<<len<<std::endl; 00222 } 00223 //zhangy add 00224 }//end while 00225 }//end rc<>rhel 00226 } else if ((rmin>gvin) && (rmax<gvout)) { 00227 // case D (contained) 00228 if (rc > rhel) { 00229 double alp = asin(rhel/rc); 00230 phim = pc - alp; 00231 phip = pc + alp; 00232 } 00233 // if rc<rhel and case D, it's an orbiter. Use max phim, phip, so do nothing 00234 } 00235 phim -= phiex; 00236 phip += phiex; 00237 00238 if(5 == debug) { 00239 std::cout<<"add hits phim "<<phim <<" phip "<<phip<< std::endl; 00240 } 00241 00242 // now try to add hits 00243 kkl = 0; 00244 if(5 == debug) std::cout<<"===== try to add hits:"<< std::endl; 00245 while(hhhh[kkl]) { 00246 MdcxHit* temphit = hhhh[kkl++]; 00247 //yzhang add 2011-10-11 00248 if(temphel->Doca_Len()<0) { 00249 if(5 == debug) std::cout<< " len<0 " << temphel->Doca_Len()<< " continue"<<std::endl; 00250 continue; 00251 } 00252 //zhangy 00253 double factor=1.; 00254 //if ((0 != temphit->type()) && (temphit->Layer()>7)) factor = MdcxParameters::addHitFactor;//yzhang FIXME 2009-11-03 00255 if(5 == debug) { 00256 temphit->print(std::cout); 00257 std::cout<< " pull "<<temphit->pull(*temphel) 00258 << " nsig "<<factor * MdcxParameters::nSigAddHitTrk<< " len "<<temphel->Doca_Len() <<std::endl; 00259 } 00260 //yzhang 2009-10-21 16:55:25 ///FIXME 00261 if((!temphit->GetUsedOnHel())&&(fabs(temphit->pull(*temphel))< factor * MdcxParameters::nSigAddHitTrk) ) { 00262 00263 //if( (!temphit->GetUsedOnHel()) && (fabs(temphit->d())<1.2) ) //delete 00264 double pw = temphit->pw(); 00265 if((phip-pw) > m_2pi) pw += m_2pi; 00266 if((phim-pw) < -m_2pi) pw -= m_2pi; 00267 if(5 == debug) { 00268 std::cout<< "phi wire "<<pw << " phim "<<phim<< " phip "<<phip<<" len "<<temphel->Doca_Len(); 00269 } 00270 if ( (pw>phim) && (pw<phip) ) { 00271 if (5 == debug){ 00272 std::cout<< " used "<<temphit->GetUsedOnHel() 00273 <<" pull " << fabs(temphit->pull(*temphel)) 00274 <<" <? nSig " << MdcxParameters::nSigAddHitTrk 00275 << std::endl; 00276 } 00277 temphit->SetConstErr(0); 00278 double pull = temphit->pull( *temphel ); 00279 ncalc++; 00280 sumpull += fabs(pull); 00281 //cout << "MdcxAddHits: hit " << kkl-1 << " trk " << whichtrack << " pull " << pull << endl; 00282 if(g_addHitCut) g_addHitCut->fill(fabs(pull)); 00283 int layer = temphit->Layer(); 00284 if(g_addHitCut2d) g_addHitCut2d->fill(layer,fabs(pull)); 00285 if(5 == debug) { 00286 std::cout<< " pull "<<pull 00287 << " addcut "<< MdcxParameters::nSigAddHitTrk 00288 << "* factor:"<< factor <<"="<<factor * MdcxParameters::nSigAddHitTrk 00289 << " len "<< temphel->Doca_Len() 00290 << " >? lmax "<< lmax 00291 << " >? maxTkLen "<< MdcxParameters::maxTrkLength 00292 << std::endl; 00293 } 00294 // yzhang 2009-10-21 22:55:26 00295 if(fabs(pull) < factor * MdcxParameters::nSigAddHitTrk) { 00296 //if(fabs(pull) < addcut) 00297 double len = temphel->Doca_Len(); 00298 //{ hhhh[kkl-1]->print(std::cout); 00299 // cout << " trk " << whichtrack << " pull " << pull 00300 // << " d() " << temphit->d() << " len " << len << endl; } 00301 //int theflag = temphel->GetTurnFlag(); //zoujh ??? 00302 if (((len>0.0)&&(len<lmax))||((len<0.0)&&(len>-MdcxParameters::maxTrkLength))) { 00303 //templist.append(temphit); 00304 listoasses.append(temphit); 00305 temphit->SetUsedOnHel(1); 00306 nadded++; 00307 if(debug>2){ 00308 temphit->print(std::cout); 00309 std::cout<<"Added "<< std::endl; 00310 } 00311 } 00312 } 00313 temphit->SetConstErr(1); 00314 }//end phim,phip cuts 00315 } 00316 } 00317 //listoasses = templist; 00318 return listoasses; 00319 }
float MdcxAddHits::addcut [protected] |
HepAList<MdcxHit> MdcxAddHits::hhhh [protected] |
int MdcxAddHits::kkk [protected] |
int MdcxAddHits::kkl [protected] |
HepAList<MdcxHit> MdcxAddHits::listoasses [protected] |
int MdcxAddHits::nadded [protected] |
int MdcxAddHits::ncalc [protected] |
int MdcxAddHits::nhits [protected] |
int MdcxAddHits::ntracks [protected] |
double MdcxAddHits::sumpull [protected] |
int MdcxAddHits::thot [protected] |
HepAList<MdcxHel> MdcxAddHits::tttt [protected] |
int MdcxAddHits::tuot [protected] |
Definition at line 54 of file MdcxAddHits.h.