GFsegment.cpp Source File

00001 
00002 
00003 //--------------------------------------------------------
00004 //
00005 //   GFsegment 
00006 //
00007 //--------------------------------------------------------
00008 #include "TkrRecon/GFsegment.h"
00009 #include "TkrRecon/GFparticle.h"
00010 
00011 //Add 9/12/01 to put in tower info
00012 #include "idents/TowerId.h"
00013 
00014 //#########################################################
00015 GFsegment::GFsegment(const GFtrack* _GFtrack)
00016 : _mGFtrack(_GFtrack)
00017 //#########################################################
00018 {
00019     m_axis = _GFtrack->getAxis(); 
00020     clear();
00021 }
00022 //#########################################################
00023 KalPlane GFsegment::getKPlane() const
00024 //#########################################################
00025 {   
00026     if (m_nextKplane.getIDPlane() < 0) {
00027         // std::cout << " error GFsegment::getKPlane " << "\n";
00028     }
00029     return m_nextKplane;
00030 }
00031 
00032 //########################################################
00033 void GFsegment::clear()
00034 //########################################################
00035 {       
00036     m_indexhit = -1;
00037     m_statusHit = GFbase::EMPTY;
00038     m_nextKplane.clear();
00039     KalTrack::clear();
00040 }
00041 //########################################################
00042 void GFsegment::next(int jplane)
00043 //########################################################
00044 {       
00045     clear();
00046     m_nextKplane.clear();
00047     
00048     KalPlane oriKplane = _mGFtrack->lastKPlane();
00049 
00050     int trkSize        = _mGFtrack->kplanelist.size();
00051 
00052     if (trkSize > 0)
00053     {
00054         if (trkSize > 1)
00055         {
00056             KalPlane prevPlane = _mGFtrack->previousKPlane();
00057             kplanelist.push_back(prevPlane);
00058         }
00059 
00060         kplanelist.push_back(oriKplane);
00061     }
00062         
00063     doit(oriKplane,jplane,KalHit::FIT);
00064     
00065 }
00066 //########################################################
00067 void GFsegment::previous(int jplane)
00068 //########################################################
00069 {       
00070     clear();
00071     m_nextKplane.clear();
00072     
00073     KalPlane oriKplane = _mGFtrack->firstKPlane();
00074 
00075     if (oriKplane.getIDPlane() != _mGFtrack->originalKPlane().getIDPlane()) 
00076         kplanelist.push_back(oriKplane);
00077     
00078     doit(oriKplane,jplane,KalHit::SMOOTH);
00079     
00080 }
00081 
00082 //#########################################################################
00083 void GFsegment::best(int klayer)
00084 //#########################################################################
00085 {
00086     int indexhit1;
00087     std::vector<int> indexhit1list;
00088     indexhit1list.clear();
00089     
00090     double sigmac = _mGFtrack->sigmaCut();
00091     double best_chiSqhit = GFcontrol::minSegmentHits*sigmac*sigmac;
00092     GFsegment best_GFsegment(_mGFtrack);
00093     //unused:   int indexhit = -1;
00094     
00095     bool found      = false;
00096         bool done       = false;
00097     bool attemptFix  = true;
00098         
00099         while (!done) 
00100     {
00101             next(klayer);
00102             
00103             double chiSq = 1e6;
00104 
00105             if (!accept()) 
00106         {
00107             if (attemptFix && numDataPoints() == 2) 
00108             {
00109                 int badHitIndex = kplanelist.back().getIDHit();
00110                         indexhit1list.push_back(badHitIndex);
00111                 GFtutor::_DATA->flagHit(m_axis,badHitIndex);
00112                 attemptFix = false;
00113             }
00114             else done = true;
00115             } 
00116         else 
00117         {
00118                     indexhit1 = m_indexhit;
00119                     chiSq     = chiGFSq();
00120 
00121                     indexhit1list.push_back(indexhit1);
00122 
00123             GFtutor::_DATA->flagHit(m_axis,indexhit1);
00124             }
00125 
00126             if (accept() && chiSq < best_chiSqhit) 
00127         {
00128                     found            = true;
00129                     best_chiSqhit    = chiSq;
00130                     KalHit hitsmooth = getKPlane(klayer).getHit(KalHit::SMOOTH);
00131                     KalHit hitfit    = m_nextKplane.getHit(KalHit::FIT);
00132                     KalHit hitNew(KalHit::FIT,hitsmooth.getPar(),hitfit.getCov());
00133                     m_nextKplane.setHit(hitNew);
00134                     best_GFsegment   = *this;
00135             }
00136         }
00137         
00138         if (!found) clear();
00139         else *this = best_GFsegment;  
00140         
00141         for (unsigned int i=0; i<indexhit1list.size(); i++) 
00142             GFtutor::_DATA->unflagHit(m_axis, indexhit1list[i]);
00143         
00144 }
00145 
00146 //#########################################################
00147 bool GFsegment::accept() const
00148 //#########################################################
00149 {
00150     bool accept = true;
00151     
00152     if (numDataPoints() < 3) {
00153         accept = false;
00154     }
00155     if (!accept) return accept;
00156     
00157     if (numGaps() > GFcontrol::maxGapsSegment) {
00158         accept = false;
00159     }
00160     if (!accept) return accept;
00161     
00162     if (chiGFSq() > GFcontrol::maxChiSqSegment) {
00163         accept = false;
00164     }
00165     if (!accept) return accept;
00166     
00167     return accept;
00168 }
00169 
00170 //########################################################
00171 void GFsegment::flagUsedHits(int jplane)
00172 //########################################################
00173 {       
00174     unsigned i = 0;
00175     for ( ; i < kplanelist.size(); i++) {
00176         int kplane = kplanelist[i].getIDPlane();
00177         int indexhit = kplanelist[i].getIDHit();
00178         if (kplane >= jplane) {
00179             GFtutor::_DATA->flagHit(m_axis,indexhit);
00180         }
00181     }
00182     
00183 }
00184 
00185 //########################################################
00186 void GFsegment::unFlagUsedHits(int jplane)
00187 //########################################################
00188 {       
00189     for (int i = 0 ; i < kplanelist.size(); i++) 
00190     {
00191             int kplane = kplanelist[i].getIDPlane();
00192             int indexhit = kplanelist[i].getIDHit();
00193         if (kplane >= jplane) 
00194         {
00195             GFtutor::_DATA->unflagHit(m_axis,indexhit);
00196         }
00197     }
00198     
00199 }
00200 //########################################################
00201 void GFsegment::unFlagAllHits()
00202 //########################################################
00203 {       
00204     for (int i = 0; i < kplanelist.size(); i++) 
00205     {
00206             //unused:           int kplane = kplanelist[i].getIDPlane();
00207             int indexhit = kplanelist[i].getIDHit();
00208             GFtutor::_DATA->unflagHit(m_axis,indexhit);
00209     }
00210     
00211 }
00212 
00213 //########################################################
00214 double GFsegment::chiGFSq() const
00215 //########################################################
00216 {       
00217     double chiGF = 1e6;
00218     if (numDataPoints() < GFcontrol::minSegmentHits) return chiGF;
00219     
00220     double sigmaC = _mGFtrack->sigmaCut();
00221     chiGF =numGaps()*sigmaC*sigmaC/(1.*GFcontrol::minSegmentHits);
00222     chiGF += chiSquareSmooth();
00223     return chiGF;
00224 }
00225 //--------------------------------------------------------
00226 // GFsegment - private
00227 //--------------------------------------------------------
00228 
00229 //#########################################################
00230 KalPlane GFsegment::getKPlane(int kplane) const
00231 //#########################################################
00232 {
00233     
00234     KalPlane Kplane;
00235     
00236     //unused:   bool done = false;
00237     unsigned i = 0;
00238     for (; i< kplanelist.size(); i++) {
00239         if (kplanelist[i].getIDPlane() == kplane ) {
00240             Kplane =  kplanelist[i];
00241             break;
00242         }
00243     }
00244     
00245     return Kplane;
00246 }
00247 //#########################################################
00248 KalPlane GFsegment::followingKPlane(int kplane) const
00249 //#########################################################
00250 {
00251     
00252     KalPlane Kplane;
00253     
00254     //unused:   bool done = false;
00255     unsigned i = 0;
00256     for (; i< kplanelist.size(); i++) {
00257         if (kplanelist[i].getIDPlane() > kplane ) {
00258             Kplane =  kplanelist[i];
00259             break;
00260         }
00261     }
00262     
00263     return Kplane;
00264     
00265 }
00266 //########################################################
00267 void GFsegment::doit(KalPlane& oriKplane, int jplane, KalHit::TYPE type)
00268 //########################################################
00269 {
00270     int iplane = oriKplane.getIDPlane();
00271     
00272     int up = ( iplane < jplane ? +1 : -1);
00273 
00274     //Loop over planes trying to form a "segment" with three hits
00275     int kplane = jplane;
00276     int lstgaps = 0;
00277     for (kplane = jplane ; -1 < kplane && kplane < GFtutor::numPlanes(); kplane+=up) 
00278     {
00279                 KalPlane prevKplane;
00280                 KalPlane nextKplane;
00281 
00282                 if (kplanelist.size() == 0) prevKplane = oriKplane; 
00283                 else                        prevKplane = kplanelist.back();
00284         
00285                 if (kplane != jplane) type = KalHit::FIT;
00286 
00287                 GFbase::StatusHit statushit = nextKPlane(prevKplane, kplane, nextKplane, type);
00288         
00289         //Found a candidate hit, include info here
00290         if (statushit == GFbase::FOUND) 
00291         { 
00292                         lstgaps = 0;
00293             
00294                         kplanelist.push_back(nextKplane);
00295                         if (kplanelist.size() >= 2) 
00296             {
00297                 filterStep(kplanelist.size()-2);
00298             }
00299                         else 
00300             {
00301                                 KalHit hitpred = nextKplane.getHit(KalHit::PRED);
00302                                 KalHit hitmeas = nextKplane.getHit(KalHit::MEAS);
00303                                 KalPar p(hitmeas.getPar().getPosition(),hitpred.getPar().getSlope());
00304                                 KalHit hitfit(KalHit::FIT,p,hitpred.getCov());
00305                                 kplanelist[0].setHit(hitfit);
00306                         }
00307                 } 
00308         //No hit found, determine what to do
00309         else
00310         {
00311             int numNew = numDataPoints() - 1;
00312             int numOld = _mGFtrack->kplanelist.size();
00313 
00314             if (numNew + numOld < 3) break;
00315 
00316             lstgaps++;
00317         }
00318         
00319                 if (kplane == jplane) 
00320         {
00321                         m_statusHit = statushit;
00322 
00323                         if (statushit != GFbase::FOUND) 
00324             {
00325                 break;
00326             }
00327 
00328                         m_nextKplane = kplanelist.back();  // save the plane - the PRED changes afther doFit()
00329                         m_indexhit       = m_nextKplane.getIDHit();
00330                 }
00331         
00332                 if (lstgaps == GFcontrol::maxConsecutiveGaps) 
00333         {
00334             break;
00335         }
00336                 if (numDataPoints() == GFcontrol::minSegmentHits) 
00337         {
00338             break;
00339         }
00340     }
00341     
00342     doFit();
00343 }
00344 //########################################################
00345 GFbase::StatusHit GFsegment::nextKPlane(const KalPlane& previousKplane, 
00346                                         int kplane, KalPlane& nextKplane,
00347                                         KalHit::TYPE type) const
00348 //########################################################
00349 {
00350     nextKplane = projectedKPlane(previousKplane, kplane, type);
00351     
00352     int indexhit = -1;
00353     double radius = 0.;
00354     GFbase::StatusHit statushit; 
00355 
00356     double sigma = sigmaFoundHit(previousKplane, nextKplane, indexhit, radius);
00357 
00358     if (sigma < _mGFtrack->m_sigmaCut) 
00359     {
00360           statushit = GFbase::FOUND;
00361           incorporateFoundHit(nextKplane, indexhit);
00362     }
00363     else statushit = GFbase::EMPTY;
00364     
00365     return statushit;
00366 }
00367 
00368 //#########################################################################
00369 KalPlane GFsegment::projectedKPlane(KalPlane prevKplane, int klayer, 
00370                                     KalHit::TYPE type) const
00371 //########################################################################
00372 {
00373     // The z end position of the next klayer plane
00374     // KalHit::TYPE type = KalHit::FIT;
00375     double zEnd = GFsegment::getZklayer(m_axis, klayer);
00376     
00377     KalHit predhit = prevKplane.predicted(type, zEnd, klayer);
00378     KalPlane projectedKplane(prevKplane.getIDHit(),klayer,prevKplane.getEnergy(), zEnd, 
00379         prevKplane.getOrthPar(), predhit);
00380     
00381     return projectedKplane;
00382 }
00383 //#########################################################################
00384 void GFsegment::incorporateFoundHit(KalPlane& nextKplane, int indexhit) const
00385 //#########################################################################
00386 {                       
00387     Point nearHit = GFtutor::_DATA->position(m_axis, indexhit);
00388     double x0 =0.;
00389     double z0 = nearHit.z();
00390     x0 = (_mGFtrack->m_axis == SiCluster::X ? nearHit.x() : nearHit.y());
00391     KalPar measpar(x0,0.);
00392     
00393     double  sigma   = GFtutor::siResolution();
00394     double  size    = GFtutor::_DATA->size(m_axis,indexhit);
00395     //int     towerId = GFtutor::_DATA->getHit(indexhit)->tower();
00396 
00397     idents::TowerId tower   = idents::TowerId(GFtutor::_DATA->getHit(indexhit)->tower());
00398     int     towerId = 10*(tower.iy()+1) + (tower.ix()+1);
00399 
00400         double factor = 1.;
00401     //if (GFcontrol::sigmaCluster) factor = 1./size; // this must be the correct one But
00402     if (GFcontrol::sigmaCluster) factor = size*size; 
00403         else factor = size*size;
00404 
00405         KalMatrix meascov(sigma*sigma*factor,0.,0.); 
00406     
00407     KalHit meashit(KalHit::MEAS, measpar, meascov);
00408     
00409     nextKplane.setIDHit(indexhit,towerId);
00410     nextKplane.setZPlane(z0);
00411     nextKplane.setHit(meashit);
00412 }
00413 //-------------------------------------------------------------
00414 //   GFsegment -> Finding the Hit
00415 //   9/18/01 Tracy Usher 
00416 //   Have tried to patch this routine a bit. Try to add a check to make
00417 //   sure hit is no further away than a neighboring tower... 
00418 //-------------------------------------------------------------
00419 //#########################################################################
00420 double GFsegment::sigmaFoundHit(const KalPlane& previousKplane, const KalPlane& nextKplane,
00421                                 int& indexhit, double& radiushit) const
00422                                 //#########################################################################
00423 {
00424     indexhit = -1;
00425     radiushit = 1e6;
00426     double sigmahit = 1e6;
00427 
00428     //Set up for current tower and tower tolerance
00429     double TowerWidth = 4.4 * GFtutor::trayWidth();
00430     int    oldIndex   = previousKplane.getIDHit();
00431     Point  oldCoords(0.,0.,0.);
00432 
00433     if (oldIndex > -1)
00434     {
00435         oldCoords  = GFtutor::_DATA->position(_mGFtrack->m_axis, oldIndex);
00436         TowerWidth = 1.1 * GFtutor::trayWidth();
00437     }
00438 
00439     
00440     double MAX_RADIUS = GFtutor::trayWidth()/2.;
00441     //double ERROR_ZPLANE= GFtutor::siThickness(); 
00442     
00443     KalHit hitp = nextKplane.getHit(KalHit::PRED);
00444 
00445     //Use predicted position to estimate the search region
00446     double xError = sqrt(hitp.getCov().getsiga());
00447     double outRadius = _mGFtrack->m_sigmaCut*xError;
00448 
00449     //The below commented out since I don't understand the zError part
00450     //I think the idea is to increase the search region for large angle 
00451     //tracks, but I am not sure why you would want to do that...
00452     //double zError=ERROR_ZPLANE*hitp.getPar().getSlope();
00453     //double rError=sqrt(xError*xError+zError*zError);
00454     //outRadius=3.*_mGFtrack->m_sigmaCut*rError;
00455     
00456     //Keep outside radius from getting too big
00457     if (outRadius > MAX_RADIUS ) outRadius = MAX_RADIUS;                
00458     
00459     double x0=0;
00460     double y0=0;
00461     if (_mGFtrack->m_axis==SiCluster::X)
00462     {
00463         x0 = hitp.getPar().getPosition();
00464         y0 = previousKplane.getOrthPar().getPosition();
00465         //y0 = oldCoords.y();
00466     }
00467     else
00468     {
00469         //x0 = oldCoords.x();
00470         x0 = previousKplane.getOrthPar().getPosition();
00471         y0 = hitp.getPar().getPosition();
00472     }
00473 
00474     double z0=nextKplane.getZPlane();
00475 
00476     Point centerX(x0,y0,z0);
00477     Point nearHit(0.,0.,z0);
00478     
00479     double inerRadius = -1.;
00480     int klayer = nextKplane.getIDPlane();
00481     double slope = hitp.getPar().getSlope();
00482     
00483     // Must be inside Glast
00484         bool done = false;
00485         while (!done) {
00486             nearHit = GFtutor::_DATA->nearestHitOutside(m_axis, klayer, inerRadius, centerX, indexhit);
00487             done = foundHit(indexhit, inerRadius, outRadius, TowerWidth, centerX, nearHit, slope);
00488         }
00489     
00490     if (indexhit >= 0) 
00491     {
00492             double x_hit        = (_mGFtrack->m_axis == SiCluster::X? nearHit.x() : nearHit.y());
00493             double x_center = (_mGFtrack->m_axis == SiCluster::X? centerX.x() : centerX.y());
00494             radiushit = fabs(x_hit-x_center);
00495             if (xError > 0.) sigmahit= radiushit/xError;
00496             //if (rError > 0.) sigmahit= radiushit/rError;
00497     }
00498     
00499     return sigmahit;
00500 }
00501 
00502 //#########################################################################
00503 bool GFsegment::foundHit(int& indexhit, double& inerRadius, double outRadius, double twrRadius,
00504                          const Point& centerX, const Point& nearHit, double slope) const
00505                          //#########################################################################
00506 {
00507     bool done = true;
00508     
00509     if (indexhit < 0) return done;
00510 
00511     double hitDelta = fabs(nearHit.x() - centerX.x());
00512     double twrDelta = fabs(nearHit.y() - centerX.y());
00513 
00514     if (_mGFtrack->m_axis == SiCluster::Y)
00515     {
00516         hitDelta = fabs(nearHit.y() - centerX.y());
00517         twrDelta = fabs(nearHit.x() - centerX.x());
00518     }
00519     
00520     //if (hitDelta < outRadius && twrDelta < twrRadius)
00521     if (hitDelta < outRadius && twrDelta < twrRadius)
00522     {
00523           if (GFtutor::_DATA->hitFlagged(_mGFtrack->m_axis, indexhit)) done = false;
00524     } 
00525     else indexhit = -1; // outside region 
00526 
00527     //double deltaX = (_mGFtrack->m_axis == SiCluster::X? fabs(nearHit.x()-centerX.x()):
00528     //fabs(nearHit.y()-centerX.y()));
00529     
00530     //if (deltaX < outRadius) {
00531         //if (GFtutor::_DATA->hitFlagged(_mGFtrack->m_axis, indexhit)) done = false;
00532     //} else indexhit = -1; // outside region 
00533     
00534     // this condition is necessary for the large angle pattern recognition
00535     if (indexhit > 0 ) 
00536     {
00537           if (GFtutor::okClusterSize(m_axis,indexhit,slope) == 0) done = false;
00538     }
00539     
00540     if (done == false) 
00541     {
00542           indexhit = -1;
00543           //inerRadius = deltaX + 0.1*GFtutor::siResolution();
00544           inerRadius = hitDelta + 0.1*GFtutor::siResolution();
00545     }
00546     
00547     return done;
00548 }
00549 //--------------------------------------------------------
00550 // GFsegment - Utility functions
00551 //--------------------------------------------------------
00552 //########################################################
00553 double GFsegment::getZklayer(enum SiCluster::view axis, int klayer) const
00554 //########################################################
00555 {
00556     // Get the zEnd of the next plane   
00557     double zEnd =0;
00558     
00559     Point oneHit = GFtutor::_DATA->meanHit(axis, klayer);
00560     zEnd = oneHit.z();
00561     
00562     if (zEnd==0. && kplanelist.size()>0) zEnd = -50.; 
00563     
00564     return zEnd;
00565 }           
00566 //########################################################
00567 bool GFsegment::crack(const KalPlane& nextKplane) const
00568 //########################################################
00569 {
00570     bool crack = false;
00571     /*
00572     KalHit hitp = nextKplane.getHit(KalHit::PRED);
00573     double x0=0.;
00574     double y0=0.;
00575     if (m_axis==SiCluster::X)  x0=hitp.getPar().getPosition();
00576     else  y0=hitp.getPar().getPosition();
00577     double z0=nextKplane.getZPlane();
00578     Point centerX(x0,y0,z0);
00579     
00580     float tower_width = GFtutor::side;
00581         + 2.*Tower::wall_thickness + Glast::wall_gap;
00582     float x_twr_bnd = fmod((centerX.x() + tower_width*2.5), tower_width);
00583     x_twr_bnd = .5*tower_width - fabs(x_twr_bnd -tower_width*.5);
00584     float y_twr_bnd = fmod((centerX.y() + tower_width*2.5), tower_width);
00585     y_twr_bnd = .5*tower_width - fabs(y_twr_bnd - tower_width*.5);
00586     float twr_bnd = (m_axis == SiCluster::X) ? x_twr_bnd : y_twr_bnd;
00587     if(twr_bnd < .5) crack = true;
00588     */
00589     return crack;
00590 }