MdcxAddHits Class Reference

#include <MdcxAddHits.h>

List of all members.

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< MdcxHeltttt
HepAList< MdcxHithhhh
HepAList< MdcxHitlistoasses


Detailed Description

Definition at line 33 of file MdcxAddHits.h.


Constructor & Destructor Documentation

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]

Definition at line 83 of file MdcxAddHits.cxx.

00083                           { 
00084 }


Member Function Documentation

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 }


Member Data Documentation

float MdcxAddHits::addcut [protected]

Definition at line 51 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

HepAList<MdcxHit> MdcxAddHits::hhhh [protected]

Definition at line 57 of file MdcxAddHits.h.

Referenced by GetAssociates(), and MdcxAddHits().

int MdcxAddHits::kkk [protected]

Definition at line 54 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

int MdcxAddHits::kkl [protected]

Definition at line 54 of file MdcxAddHits.h.

Referenced by GetAssociates(), and MdcxAddHits().

HepAList<MdcxHit> MdcxAddHits::listoasses [protected]

Definition at line 58 of file MdcxAddHits.h.

Referenced by GetAssociates().

int MdcxAddHits::nadded [protected]

Definition at line 52 of file MdcxAddHits.h.

Referenced by GetAssociates().

int MdcxAddHits::ncalc [protected]

Definition at line 52 of file MdcxAddHits.h.

Referenced by GetAssociates().

int MdcxAddHits::nhits [protected]

Definition at line 55 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

int MdcxAddHits::ntracks [protected]

Definition at line 55 of file MdcxAddHits.h.

Referenced by GetAssociates(), and MdcxAddHits().

double MdcxAddHits::sumpull [protected]

Definition at line 53 of file MdcxAddHits.h.

Referenced by GetAssociates().

int MdcxAddHits::thot [protected]

Definition at line 54 of file MdcxAddHits.h.

Referenced by MdcxAddHits().

HepAList<MdcxHel> MdcxAddHits::tttt [protected]

Definition at line 56 of file MdcxAddHits.h.

Referenced by GetAssociates(), and MdcxAddHits().

int MdcxAddHits::tuot [protected]

Definition at line 54 of file MdcxAddHits.h.


Generated on Tue Nov 29 23:20:20 2016 for BOSS_7.0.2 by  doxygen 1.4.7