GlastFit.cxx Source File

00001 // $Id: GlastFit.cxx,v 1.11 2000/12/11 16:38:45 burnett Exp $
00002 //------------------------------------------------------------------------------
00003 //
00004 //     GlastFit
00005 //
00006 //             Is a Kalman Track Follower class
00007 //
00008 //             It uses a Ray direction as starting seed
00009 //             and picks Si-Clusters inside a n-sigma region of the projected track into the Si Plane
00010 //
00011 //
00012 //                              B. Atwood
00013 //                              B. Atwood, JA Hernando    Santa Cruz 2/26/99
00014 //------------------------------------------------------------------------------
00015 // The following is needed because instrument/REQ got included
00016 
00017 #include "reconstruction/GlastFit.h"
00018 
00019 #include "reconstruction/SiClusters.h"
00020 
00021 
00022 static int number_of_trays = 19;
00023 static double  mod_width        = 19.485000;
00024 static double  si_thickness     = 0.04000;
00025 static double  si_strip_pitch   = 0.020848214285714;
00026 
00027 static double tower_mod_width = 38.970000000000;
00028 static double tower_wall_thickness = 0.15000000000000;
00029 static double glast_wall_gap = 0.10000000000000;
00030 
00031 #include "geometry/Ray.h"
00032 #include <iomanip>
00033 #include <algorithm>
00034 #include <cmath>
00035 using std::min;
00036 //
00037 // -------- Update TrackerRecon
00038 //
00039 
00040 static int                 CONTROL_error              = 0;
00041 static int                 CONTROL_maxConsecutiveGaps = 2;      // max consecutive Gaps - Stop
00042 static int                 CONTROL_minSegmentHits     = 3;      // min number of hits for segment
00043 static double              CONTROL_minEnergy          = 0.03;   // min tracking energy GeV
00044 static double              CONTROL_xEne               = 0.50;   // initial sharing of the energy
00045 
00046 static double       CONTROL_iniErrorSlope = 0.34; // 20 deg
00047 static double       CONTROL_iniErrorPosition = 0.01; // 0.1 mm
00048 
00049 static double CUT_maxChiSqSegment     = 200.;    // max chi2 of the initial segment
00050 static double CUT_maxGapsSegment      = 1;     // gaps allowed in the initial segment
00051 static double CUT_maxSigmasSharedHits = 1./4.;     // allowed share it it whitin nsigmas
00052 static double CUT_maxChiSq            = 1e6;
00053 static double CUT_minQuality          = -1e6;
00054 
00055 static double CUT_sigmaVeto = 8.;
00056 
00057 static double CONTROL_maxSigmaCut = 8.;
00058 
00059 //unused: static double CUT_maxSigmaSharing    = 4;
00060 //unused: static double CUT_numSigmaCutPairTrack   = 8.;
00061 //unused: static double CUT_numSigmaCutBestTrack   = 8.;
00062 
00063 //unused: static bool CONTROL_simpleStep = false;
00064 //unused: static bool CONTROL_selfishStep = true;
00065 static bool CONTROL_addTracksChi2 = true;
00066 
00067 //unused: static int CONTROL_writeOut = 2;
00068 
00069 SiClusters* GFtutor::_DATA = 0;// Sets by TrackerRecon - made accesible to all GF-objects
00070 
00071 bool GFtutor::CUT_veto = true;
00072 bool GFtutor::CONTROL_connectGFpair = true;
00073 //------------------------------------------------------------------------------
00074 //
00075 //     GFtrack
00076 //
00077 //             Is a Kalman Track Follower class
00078 //
00079 //             It uses a Ray direction as starting seed
00080 //             and picks Si-Clusters inside a n-sigma region of the projected track into the Si Plane
00081 //
00082 //
00083 //                              B. Atwood
00084 //                              B. Atwood, JA Hernando    Santa Cruz 2/26/99
00085 //------------------------------------------------------------------------------
00086 
00087 //##########################################
00088 GFbase::GFbase(double sigmaCut, double ene, int ist, const Ray& testRay):
00089 m_sigmaCut(sigmaCut),
00090 m_iniEnergy(ene),
00091 m_iniLayer(ist)
00092 //##########################################
00093 {
00094     // control defauls
00095     m_alive = true;
00096 
00097     // input data
00098     if (m_iniEnergy < CONTROL_minEnergy) m_iniEnergy = CONTROL_minEnergy;
00099     m_inVertex = testRay.position();
00100     m_inDirection = testRay.direction();
00101 
00102 }       
00103 //########################################################
00104 void GFbase::doit()                             
00105 //########################################################
00106 {
00107     int kplane = m_iniLayer;
00108 
00109     for ( ; -1 < kplane && kplane < number_of_trays-1; kplane++){
00110 
00111         step(kplane);
00112         anastep(kplane);
00113         if (!alive()) {
00114             break;
00115         }
00116     }
00117 }
00118 
00119 //##########################################
00120 void GFdata::ini()
00121 //##########################################
00122 {
00123     m_vertex = Point(0.,0.,0.);
00124     m_direction = Vector(0.,0.,0.);
00125     m_RCenergy=0.;
00126     m_quality=-100.;
00127     m_firstLayer=-1;
00128     m_nhits = -1;
00129     m_itower =-1;
00130 }
00131 //##########################################
00132 bool GFdata::empty() const
00133 //##########################################
00134 {
00135     bool empty = false;
00136     if (m_firstLayer < 0) empty = true;
00137     return empty;
00138 }
00139 
00140 //########################################################
00141 void GFdata::writeOut(int level) const
00142 //########################################################
00143 {
00144     std::cout << " --- GFdata::writeOut --- " << "\n";
00145     std::cout << " Quality       = " << Q() << "\n";
00146     std::cout << " Vertex        = " << vertex().x() << " " <<  vertex().y() << " " << vertex().z() << "\n";
00147     std::cout << " Direction     = " << direction().x() << " " << direction().y() << " " << direction().z() << "\n";
00148     std::cout << " RCenergy      = " << RCenergy() << "\n";
00149     std::cout << " first Layer   = " << firstLayer() << "\n";
00150     std::cout << " Tower         = " << tower() << "\n";
00151     std::cout << " num Hits      = " << nhits() << "\n";
00152 }
00153 
00154 //#########################################################################
00155 Point GFdata::doVertex(const Ray& r1, const Ray& r2)
00156 //#########################################################################
00157 {
00158     Point Vertex(0.,0.,0.);
00159 
00160     double z1 = r1.position().z();
00161     double z2 = r2.position().z();
00162     double cosz1 = r1.direction().z();
00163     double cosz2 = r2.direction().z();
00164     if (cosz2 == 0. || cosz1 == 0.) return Vertex;
00165 
00166     // comon position z2ref = z1
00167     Vector Vdir = r2.direction();
00168     Point Porigin = r2.position((z1-z2)/cosz2);
00169 
00170     Ray r2ref(Porigin,Vdir);
00171 
00172     double deltaX = r1.position().x()-r2ref.position().x();
00173     double deltaY = r1.position().y()-r2ref.position().y();
00174     double deltaSlopeX = (r2ref.direction().x()/cosz2) - (r1.direction().x()/cosz1);
00175     double deltaSlopeY = (r2ref.direction().y()/cosz2) - (r1.direction().y()/cosz1);
00176 
00177     bool xview = true;
00178     bool yview = true;
00179     double deltaZX = 0.;
00180     double deltaZY = 0.;
00181     if (deltaSlopeX == 0.) xview = false;
00182     if (deltaSlopeY == 0.) yview = false;
00183     if (xview) deltaZX = deltaX/deltaSlopeX;
00184     if (yview) deltaZY = deltaY/deltaSlopeY;
00185 
00186     if (xview && !yview) Vertex = r1.position(deltaZX/cosz1);
00187     else if (!xview && yview) Vertex = r1.position(deltaZY/cosz1);
00188     else if (xview && yview) {
00189         Vertex  = (z1 >= z2? r1.position(0.) : r2.position(0.));
00190         // We assume r1 = X ray and r2 = Y ray;
00191         if (Vertex.x() == 0.) {
00192             double x0 = r1.position((z2-z1)/cosz1).x();
00193             Vertex = Point( x0, r2.position().y(),r2.position().z());
00194         }
00195         if (Vertex.y() == 0.) {
00196             double y0 = r2.position((z1-z2)/cosz2).y();
00197             Vertex = Point(r1.position().x(), y0 , r1.position().z());
00198         }
00199     }
00200     return Vertex;
00201 }
00202 
00203 //#########################################################################
00204 Vector GFdata::doDirection(const Vector& xdir, const Vector& ydir)
00205 //#########################################################################
00206 {
00207     // we assume xdir goes in the X direction and ydir goes in the Y direction
00208     double slopeX = xdir.x()/xdir.z();
00209     double slopeY = ydir.y()/ydir.z();
00210 
00211     double factor = -1;
00212     Vector dir = Vector(factor*slopeX,factor*slopeY,factor).unit();
00213 
00214     return dir;
00215 }
00216 
00217 //--------------------------------------------------------
00218 //
00219 //   GFsegment
00220 //
00221 //--------------------------------------------------------
00222 //#########################################################
00223 GFsegment::GFsegment(const GFtrack* _GFtrack)
00224 : _mGFtrack(_GFtrack)
00225 //#########################################################
00226 {
00227     m_axis = _GFtrack->getAxis();
00228     clear();
00229 }
00230 //#########################################################
00231 KalPlane GFsegment::getKPlane() const
00232 //#########################################################
00233 {
00234     if (m_nextKplane.getIDPlane() < 0) {
00235         // std::cout << " error GFsegment::getKPlane " << "\n";
00236     }
00237     return m_nextKplane;
00238 }
00239 
00240 //########################################################
00241 void GFsegment::clear()
00242 //########################################################
00243 {       
00244     m_indexhit = -1;
00245     m_statusHit = GFbase::EMPTY;
00246     m_nextKplane.clear();
00247     KalTrack::clear();
00248 }
00249 //########################################################
00250 void GFsegment::next(int jplane)
00251 //########################################################
00252 {       
00253     clear();
00254     m_nextKplane.clear();
00255 
00256     KalPlane oriKplane = _mGFtrack->lastKPlane();
00257 
00258     doit(oriKplane,jplane,KalHit::FIT);
00259 
00260 }
00261 //########################################################
00262 void GFsegment::previous(int jplane)
00263 //########################################################
00264 {       
00265     clear();
00266     m_nextKplane.clear();
00267 
00268     KalPlane oriKplane = _mGFtrack->firstKPlane();
00269 
00270     doit(oriKplane,jplane,KalHit::SMOOTH);
00271 
00272 }
00273 
00274 //#########################################################################
00275 void GFsegment::best(int klayer)
00276 //#########################################################################
00277 {
00278     int indexhit1;
00279     std::vector<int> indexhit1list;
00280     indexhit1list.clear();
00281 
00282     double sigmac = _mGFtrack->sigmaCut();
00283     double best_chiSqhit = CONTROL_minSegmentHits*sigmac*sigmac;
00284     GFsegment best_GFsegment(_mGFtrack);
00285     //unused:   int indexhit = -1;
00286 
00287     bool found = false;
00288         bool done = false;
00289         
00290         while (!done) {
00291 
00292             next(klayer);
00293         
00294             double chiSq = 1e6;
00295             if (!accept()) {
00296                 done = true;
00297             } else {
00298                 indexhit1 = m_indexhit;
00299                 chiSq = chiGFSq();
00300                 indexhit1list.push_back(indexhit1);
00301 
00302                 GFtutor::_DATA->flagHit(m_axis,indexhit1);
00303             }
00304 
00305             if (accept() && chiSq < best_chiSqhit) {
00306                 
00307                 found = true;
00308                 best_chiSqhit = chiSq;
00309                 KalHit hitsmooth = getKPlane(klayer).getHit(KalHit::SMOOTH);
00310                 KalHit hitfit = m_nextKplane.getHit(KalHit::FIT);
00311                 KalHit hitNew(KalHit::FIT,hitsmooth.getPar(),hitfit.getCov());
00312                 m_nextKplane.setHit(hitNew);
00313                 best_GFsegment = *this;
00314             }
00315         }
00316         
00317         if (!found) clear();
00318         else *this = best_GFsegment;
00319         
00320         for (unsigned int i=0; i<indexhit1list.size(); i++)
00321             GFtutor::_DATA->unflagHit(m_axis, indexhit1list[i]);
00322         
00323             /*
00324             bool found = false;
00325             bool done = false;
00326             while (!done) {
00327         
00328               next(klayer);
00329         
00330                 double chiSqhit = 1e6;
00331                 
00332                   if (status() == GFbase::EMPTY) {
00333                   done = true;
00334                   continue;
00335                   }
00336                 
00337                     double next = true;
00338                     if (numDataPoints() < CONTROL_minSegmentHits) next = false;
00339                     else {
00340                     if (followingKPlane(klayer).getIDPlane() != klayer+1) next = false;
00341                     }
00342                 
00343                       if (!next) {
00344                       indexhit1 = m_indexhit;
00345                       GFtutor::_DATA->flagHit(m_axis,indexhit1);
00346                       indexhit1list.push_back(indexhit1);
00347                 
00348                         for (unsigned int i=0; i<indexhit2list.size(); i++)
00349                         GFtutor::_DATA->unflagHit(m_axis, indexhit2list[i]);
00350                         
00351                           indexhit2list.clear();
00352                           } else {
00353                           indexhit1 = m_indexhit;
00354                         
00355                             chiSqhit = chiGFSq();
00356                             indexhit2 = followingKPlane(klayer).getIDHit();                     
00357                             indexhit2list.push_back(indexhit2);
00358                             GFtutor::_DATA->flagHit(m_axis,indexhit2);
00359                             }
00360                         
00361                               if (accept() && chiSqhit < best_chiSqhit ) {
00362                               found = true;
00363                               best_chiSqhit = chiSqhit;
00364                               KalHit hitsmooth = getKPlane(klayer).getHit(KalHit::SMOOTH);
00365                               KalHit hitfit = m_nextKplane.getHit(KalHit::FIT);
00366                               KalHit hitNew(KalHit::FIT,hitsmooth.getPar(),hitfit.getCov());
00367                               m_nextKplane.setHit(hitNew);
00368                               best_GFsegment = *this;
00369                               }
00370                               }
00371                         
00372                                 if (!found) clear();
00373                                 else *this = best_GFsegment;
00374                                 
00375                                   for (unsigned int i=0; i<indexhit1list.size(); i++)
00376                                   GFtutor::_DATA->unflagHit(m_axis, indexhit1list[i]);
00377                                   for (i=0; i<indexhit2list.size(); i++)
00378         GFtutor::_DATA->unflagHit(m_axis, indexhit2list[i]); */
00379 }
00380 
00381 //#########################################################
00382 bool GFsegment::accept() const
00383 //#########################################################
00384 {
00385     bool accept = true;
00386 
00387     if (numDataPoints() < 3) {
00388         accept = false;
00389     }
00390     if (!accept) return accept;
00391 
00392     if (numGaps() > CUT_maxGapsSegment) {
00393         accept = false;
00394     }
00395     if (!accept) return accept;
00396 
00397     if (chiGFSq() > CUT_maxChiSqSegment) {
00398         accept = false;
00399     }
00400     if (!accept) return accept;
00401 
00402     return accept;
00403 }
00404 
00405 //########################################################
00406 void GFsegment::flagUsedHits(int jplane)
00407 //########################################################
00408 {       
00409     unsigned i = 0;
00410     for ( ; i < kplanelist.size(); i++) {
00411         int kplane = kplanelist[i].getIDPlane();
00412         int indexhit = kplanelist[i].getIDHit();
00413         if (kplane >= jplane) {
00414             GFtutor::_DATA->flagHit(m_axis,indexhit);
00415         }
00416     }
00417 
00418 }
00419 //########################################################
00420 void GFsegment::unFlagAllHits()
00421 //########################################################
00422 {       
00423     unsigned i = 0;
00424     for ( ; i < kplanelist.size(); i++) {
00425         //unused:               int kplane = kplanelist[i].getIDPlane();
00426         int indexhit = kplanelist[i].getIDHit();
00427         GFtutor::_DATA->unflagHit(m_axis,indexhit);
00428     }
00429 
00430 }
00431 
00432 //########################################################
00433 double GFsegment::chiGFSq() const
00434 //########################################################
00435 {       
00436     double chiGF = 1e6;
00437     if (numDataPoints() < CONTROL_minSegmentHits) return chiGF;
00438 
00439     double sigmaC = _mGFtrack->sigmaCut();
00440     chiGF =numGaps()*sigmaC*sigmaC/(1.*CONTROL_minSegmentHits);
00441     chiGF += chiSquareSmooth();
00442     return chiGF;
00443 }
00444 //--------------------------------------------------------
00445 // GFsegment - private
00446 //--------------------------------------------------------
00447 
00448 //#########################################################
00449 KalPlane GFsegment::getKPlane(int kplane) const
00450 //#########################################################
00451 {
00452 
00453     KalPlane Kplane;
00454 
00455     //unused:   bool done = false;
00456     unsigned i = 0;
00457     for (; i< kplanelist.size(); i++) {
00458         if (kplanelist[i].getIDPlane() == kplane ) {
00459             Kplane =  kplanelist[i];
00460             break;
00461         }
00462     }
00463 
00464     return Kplane;
00465 }
00466 //#########################################################
00467 KalPlane GFsegment::followingKPlane(int kplane) const
00468 //#########################################################
00469 {
00470 
00471     KalPlane Kplane;
00472 
00473     //unused:   bool done = false;
00474     unsigned i = 0;
00475     for (; i< kplanelist.size(); i++) {
00476         if (kplanelist[i].getIDPlane() > kplane ) {
00477             Kplane =  kplanelist[i];
00478             break;
00479         }
00480     }
00481 
00482     return Kplane;
00483 
00484 }
00485 //########################################################
00486 void GFsegment::doit(KalPlane& oriKplane, int jplane, KalHit::TYPE type)
00487 //########################################################
00488 {       
00489     if (oriKplane.getIDPlane() != _mGFtrack->originalKPlane().getIDPlane())
00490         kplanelist.push_back(oriKplane);
00491 
00492     int iplane = oriKplane.getIDPlane();
00493 
00494     int up = ( iplane < jplane ? +1 : -1);
00495 
00496     int kplane = jplane;
00497     int lstgaps = 0;
00498     for ( ; -1 < kplane && kplane < number_of_trays -1; kplane+=up) {
00499         
00500         KalPlane prevKplane;
00501         KalPlane nextKplane;
00502         if (kplanelist.size() == 0) prevKplane = oriKplane;
00503         else prevKplane = kplanelist.back();
00504         
00505         if (kplane != jplane) type = KalHit::FIT;
00506         GFbase::StatusHit statushit = nextKPlane(prevKplane, kplane, nextKplane, type);
00507         
00508         if (statushit == GFbase::FOUND) {
00509             lstgaps = 0;
00510 
00511             kplanelist.push_back(nextKplane);
00512             if (kplanelist.size() >= 2) filterStep(kplanelist.size()-2);
00513             else {
00514                 KalHit hitpred = nextKplane.getHit(KalHit::PRED);
00515                 KalHit hitmeas = nextKplane.getHit(KalHit::MEAS);
00516                 KalPar p(hitmeas.getPar().getPosition(),hitpred.getPar().getSlope());
00517                 KalHit hitfit(KalHit::FIT,p,hitpred.getCov());
00518                 kplanelist[0].setHit(hitfit);
00519             }
00520         } else lstgaps++;
00521         
00522         if (kplane == jplane) {
00523             m_statusHit = statushit;
00524             if (statushit != GFbase::FOUND) {
00525                 break;
00526             }
00527             m_nextKplane = kplanelist.back();  // save the plane - the PRED changes afther doFit()
00528             m_indexhit   = m_nextKplane.getIDHit();
00529         }
00530         
00531         if (lstgaps == CONTROL_maxConsecutiveGaps) {
00532             break;
00533         }
00534         if (numDataPoints() == CONTROL_minSegmentHits) {
00535             break;
00536         }
00537     }
00538 
00539     doFit();
00540 }
00541 //########################################################
00542 GFbase::StatusHit GFsegment::nextKPlane(const KalPlane& previousKplane,
00543                                         int kplane, KalPlane& nextKplane,
00544                                         KalHit::TYPE type) const
00545 //########################################################
00546 {
00547     nextKplane = projectedKPlane(previousKplane, kplane, type);
00548 
00549     int indexhit = -1;
00550     double radius = 0.;
00551     GFbase::StatusHit statushit;
00552     double sigma = sigmaFoundHit(previousKplane, nextKplane, indexhit, radius);
00553     if (sigma < _mGFtrack->m_sigmaCut) {
00554         statushit = GFbase::FOUND;
00555         incorporateFoundHit(nextKplane, indexhit);
00556     } else statushit = GFbase::EMPTY;
00557 
00558     return statushit;
00559 }
00560 
00561 //#########################################################################
00562 KalPlane GFsegment::projectedKPlane(KalPlane prevKplane, int klayer,
00563                                     KalHit::TYPE type) const
00564 //#########################################################################
00565 {
00566     // The z end position of the next klayer plane
00567     // KalHit::TYPE type = KalHit::FIT;
00568     double zEnd = GFsegment::getZklayer(m_axis, klayer);
00569 
00570     KalHit predhit = prevKplane.predicted(type, zEnd, klayer);
00571     KalPlane projectedKplane(prevKplane.getIDHit(),klayer,prevKplane.getEnergy(), zEnd,
00572         prevKplane.getOrthPar(), predhit);
00573 
00574     return projectedKplane;
00575 }
00576 //#########################################################################
00577 void GFsegment::incorporateFoundHit(KalPlane& nextKplane, int indexhit) const
00578 //#########################################################################
00579 {                       
00580     Point nearHit = GFtutor::_DATA->clusterPosition(m_axis, indexhit);
00581     double x0 =0.;
00582     double z0 = nearHit.z();
00583     x0 = (_mGFtrack->m_axis == SiData::X ? nearHit.x() : nearHit.y());
00584     KalPar measpar(x0,0.);
00585 
00586     double sigma = si_strip_pitch/sqrt(12.);
00587     double size  = GFtutor::_DATA->clusterSize(m_axis,indexhit);
00588     // KalMatrix meascov(sigma*sigma/size,0.,0.); // This must be the correct one But
00589     KalMatrix meascov(sigma*sigma*size*size,0.,0.);
00590 
00591     KalHit meashit(KalHit::MEAS, measpar, meascov);
00592 
00593     nextKplane.setIDHit(indexhit);
00594     nextKplane.setZPlane(z0);
00595     nextKplane.setHit(meashit);
00596 }
00597 //-------------------------------------------------------------
00598 //   GFsegment -> Finding the Hit
00599 //-------------------------------------------------------------
00600 //#########################################################################
00601 double GFsegment::sigmaFoundHit(const KalPlane& previousKplane, const KalPlane& nextKplane,
00602                                 int& indexhit, double& radiushit) const
00603                                 //#########################################################################
00604 {
00605     indexhit = -1;
00606     radiushit = 1e6;
00607     double sigmahit = 1e6;
00608 
00609     double DELTACHI2_CUT = _mGFtrack->m_sigmaCut;
00610     double MAX_RADIUS = mod_width/2.;
00611     double ERROR_ZPLANE = si_thickness;
00612 
00613     KalHit hitp = nextKplane.getHit(KalHit::PRED);
00614     double xError = sqrt(hitp.getCov().getsiga());
00615     double outRadius = DELTACHI2_CUT*xError;
00616     double zError=ERROR_ZPLANE*hitp.getPar().getSlope();
00617     double rError=sqrt(xError*xError+zError*zError);
00618     outRadius=3.*DELTACHI2_CUT*rError;
00619     if (outRadius > MAX_RADIUS ) outRadius = MAX_RADIUS;                
00620 
00621     double x0=0;
00622     double y0=0;
00623     if (_mGFtrack->m_axis==SiData::X)  x0=hitp.getPar().getPosition();
00624     else  y0=hitp.getPar().getPosition();
00625     double z0=nextKplane.getZPlane();
00626     Point centerX(x0,y0,z0);
00627     Point nearHit(0.,0.,z0);
00628 
00629     double inerRadius = -1.;
00630     int klayer = nextKplane.getIDPlane();
00631     double slope = hitp.getPar().getSlope();
00632 #if 0 //THB: need to do this some other way
00633 
00634     // Must be inside Glast
00635     GlastMed* _glast = GlastMed::instance();
00636     const Medium* _med0= _glast->inside(centerX);
00637     if (_med0) {
00638 #endif
00639         bool done = false;
00640         int loop_check=10000;
00641         while (!done) {
00642             nearHit = GFtutor::_DATA->nearestHitOutside(m_axis, klayer, inerRadius, centerX, &indexhit);
00643             done = foundHit(indexhit, inerRadius, outRadius, centerX, nearHit, slope);
00644                 if(--loop_check<=0){ 
00645                         std::cerr << "In GlastFit::sigmaFoundHit: infinite loop?" << std::endl; break;}
00646         }
00647 #if 0
00648     }
00649 #endif
00650     if (indexhit >= 0) {
00651         double x_hit    = (_mGFtrack->m_axis == SiData::X? nearHit.x() : nearHit.y());
00652         double x_center = (_mGFtrack->m_axis == SiData::X? centerX.x() : centerX.y());
00653         radiushit = abs(x_hit-x_center);
00654         if (rError > 0.) sigmahit= radiushit/rError;
00655     }
00656 
00657     return sigmahit;
00658 }
00659 
00660 //#########################################################################
00661 bool GFsegment::foundHit(int& indexhit, double& inerRadius, double outRadius,
00662                          const Point& centerX, const Point& nearHit, double slope) const
00663                          //#########################################################################
00664 {
00665     bool done = true;
00666 
00667     if (indexhit < 0) return done;
00668     double deltaX = (_mGFtrack->m_axis == SiData::X? std::abs(nearHit.x()-centerX.x()):
00669         std::abs(nearHit.y()-centerX.y()));
00670 
00671     if (deltaX < outRadius) {
00672         if (GFtutor::_DATA->hitFlagged(_mGFtrack->m_axis, indexhit)) done = false;
00673     } else indexhit = -1; // outside region
00674 
00675     // this condition is necessary for the large angle pattern recognition
00676     if (indexhit > 0 ) {
00677         if (GFtutor::okClusterSize(m_axis,indexhit,slope) == 0) done = false;
00678     }
00679 
00680     if (done == false) {
00681         indexhit = -1;
00682         inerRadius = deltaX + 0.1 * si_strip_pitch/sqrt(12.);
00683     }
00684 
00685     return done;
00686 }
00687 //--------------------------------------------------------
00688 // GFsegment - Utility functions
00689 //--------------------------------------------------------
00690 //########################################################
00691 double GFsegment::getZklayer(enum SiData::Axis axis, int klayer) const
00692 //########################################################
00693 {
00694     // Get the zEnd of the next plane   
00695     double zEnd =0;
00696 
00697     float rms;
00698     Point oneHit = GFtutor::_DATA->meanHit(axis, klayer, &rms);
00699     zEnd = oneHit.z();
00700 
00701     if (zEnd==0. && kplanelist.size()>0) zEnd = -50.;
00702 
00703     return zEnd;
00704 }       
00705 //########################################################
00706 bool GFsegment::crack(const KalPlane& nextKplane) const
00707 //########################################################
00708 {
00709     bool crack = false;
00710 
00711     KalHit hitp = nextKplane.getHit(KalHit::PRED);
00712     double x0=0.;
00713     double y0=0.;
00714     if (m_axis==SiData::X)  x0=hitp.getPar().getPosition();
00715     else  y0=hitp.getPar().getPosition();
00716     double z0=nextKplane.getZPlane();
00717     Point centerX(x0,y0,z0);
00718 
00719     float tower_width = mod_width
00720         + 2.*tower_wall_thickness + glast_wall_gap;
00721     float x_twr_bnd = fmod((centerX.x() + tower_width*2.5), tower_width);
00722     x_twr_bnd = .5*tower_width - fabs(x_twr_bnd -tower_width*.5);
00723     float y_twr_bnd = fmod((centerX.y() + tower_width*2.5), tower_width);
00724     y_twr_bnd = .5*tower_width - fabs(y_twr_bnd - tower_width*.5);
00725     float twr_bnd = (m_axis == SiData::X) ? x_twr_bnd : y_twr_bnd;
00726     if(twr_bnd < .5) crack = true;
00727 
00728     return crack;
00729 }
00730 
00731 //-----------------------------------------------------
00732 //
00733 //   GFtrack
00734 //
00735 //-----------------------------------------------------
00736 //######################################################
00737 GFtrack::GFtrack(SiData::Axis axis,
00738                  double sigmaCut,
00739                  double energy,
00740                  int ist,
00741                  const Ray& testRay,
00742                  bool run)
00743   : GFbase(sigmaCut,energy,ist,testRay),
00744     m_axis(axis), m_status(EMPTY),
00745     m_qbest(-1e6), m_gaps(0), m_istGaps(0), m_lstLayer(0), m_noisyHits(0),
00746     m_istNoisyHits(0)
00747 //#######################################################
00748 {
00749     ini();
00750     if (run == true) {
00751         doit();
00752         fit();
00753         if (empty()) clear();
00754     }
00755 }
00756 
00757 //##########################################
00758 void GFtrack::flagAllHits()
00759 //##########################################
00760 {
00761     for(unsigned int i=0; i<kplanelist.size(); i++) {
00762         GFtutor::_DATA->flagHit(m_axis, kplanelist[i].getIDHit());
00763     }
00764 }
00765 
00766 //##########################################
00767 void GFtrack::unFlagAllHits()
00768 //##########################################
00769 {
00770     for (unsigned i=0; i<kplanelist.size(); i++) {
00771         GFtutor::_DATA->unflagHit(m_axis, kplanelist[i].getIDHit());
00772     }
00773 }
00774 //########################################################
00775 bool GFtrack::empty() const
00776 //########################################################
00777 {
00778     bool empty = GFdata::empty();
00779     if (firstLayer() < 0) empty = true;
00780     if (numDataPoints() < CONTROL_minSegmentHits) empty = true;
00781     if (chiSquare() < 0.) empty = true;
00782     return empty;
00783 }
00784 
00785 //########################################################
00786 bool GFtrack::accept() const
00787 //########################################################
00788 {
00789     bool valid = false;
00790     if (empty()) return valid;
00791 
00792     if (chiSquare() > CUT_maxChiSq) return valid;
00793     if (Q() < CUT_minQuality) return valid;
00794     int idhit;
00795     double sigma;
00796     if (GFtutor::CUT_veto) {
00797         if (veto(idhit,sigma)) return valid;
00798     }
00799 
00800     valid = true;
00801     return valid;
00802 }
00803 //##########################################
00804 void GFtrack::clear()
00805 //##########################################
00806 {
00807 
00808     KalTrack::clear();
00809 
00810     m_lstGaps = 0;
00811     m_gaps   = 0;
00812     m_istGaps= 0;
00813     m_lstLayer = 0;
00814     m_noisyHits = 0;
00815     m_istNoisyHits = 0;
00816 
00817     m_status = EMPTY;
00818 
00819     if (_mGFsegment) _mGFsegment->clear();
00820 
00821     m_qbest = -1e6;
00822     GFdata::ini();
00823     setAlive();
00824 
00825 }
00826 //########################################################
00827 void GFtrack::writeOut(int level) const
00828 //########################################################
00829 {
00830     // kludge to avoid egcs warning messages
00831     int axis = getAxis();
00832     int status = m_status;
00833     std::cout << " --- GFtrack::writeOut --- " << "\n";
00834     std::cout << " axis           = " << axis << "\n";
00835     std::cout << " Qbest          = " << Qbest() << "\n";
00836     std::cout << " last  Layer    = " << lastLayer() << "\n";
00837     std::cout << " num Hits       = " << numDataPoints() << "\n";
00838     std::cout << " num Gaps       = " << numGaps() << "\n";
00839     std::cout << " num First Gaps = " << numFirstGaps() << "\n";
00840     std::cout << " num Noise      = " << numNoise() << "\n";
00841     std::cout << " num First Noise= " << numFirstNoise() << "\n";
00842     std::cout << " last Status    = " << status << "\n";
00843 
00844     GFdata::writeOut(level);
00845 
00846     std::cout << " --> KalTrack : " << "\n";
00847     KalTrack::writeOut(level);
00848 
00849 }
00850 
00851 //#########################################################################
00852 bool GFtrack::veto(int& idhit, double& sigma) const
00853 //#########################################################################
00854 {
00855     // WORK
00856     bool veto = false;
00857 
00858     int klayer = firstKPlane().getIDPlane() -1;
00859     if (klayer < 0) return veto;
00860     _mGFsegment->previous(klayer);
00861 
00862     if (_mGFsegment->status() == FOUND) {
00863         idhit = _mGFsegment->indexhit();
00864         sigma = _mGFsegment->getKPlane().getSigma(KalHit::PRED);
00865         if (sigma < CUT_sigmaVeto) veto = true;
00866     }
00867 
00868     return veto;
00869 }
00870 
00871 //--------------------------------------------------------
00872 //  GFtrack - Private
00873 //--------------------------------------------------------
00874 
00875 //########################################################
00876 void GFtrack::ini()
00877 //########################################################
00878 {
00879     m_status = EMPTY;
00880 
00881     m_gaps   = 0;
00882     m_istGaps= 0;
00883     m_lstGaps = 0;
00884     m_noisyHits = 0;
00885     m_istNoisyHits   = 0;
00886     m_lstLayer = 0;
00887 
00888     KalTrack::clear();
00889 
00890     _mGFsegment = 0;
00891     _mGFsegment = new GFsegment(this);
00892 
00893     m_qbest = -1e6;
00894     GFdata::ini();
00895     setAlive();
00896 
00897 }
00898 
00899 //#######################################################
00900 void GFtrack::step(int kplane)
00901 //#######################################################
00902 {
00903     if (!alive()) return;
00904 
00905     m_status = EMPTY;
00906     if (kplane > number_of_trays -1) return;
00907 
00908     if (numDataPoints() == 0) _mGFsegment->best(kplane);
00909     else {
00910         _mGFsegment->next(kplane);
00911     }
00912 
00913     m_status = _mGFsegment->status();
00914     if (m_status == FOUND) kplanelist.push_back(_mGFsegment->getKPlane());
00915 
00916 }
00917 
00918 //#########################################################
00919 void GFtrack::anastep(int kplane)
00920 //#########################################################
00921 {
00922     if (!alive()) return;
00923 
00924     contability(kplane);
00925     if (end()) {
00926         kill();
00927     }
00928 }
00929 
00930 //##########################################################
00931 void GFtrack::contability(int kplane)
00932 //##########################################################
00933 {
00934     if (!alive()) return;
00935 
00936     if (m_status != FOUND) {
00937         m_lstGaps++;
00938         if (numDataPoints()>0 && numDataPoints()<3) m_istGaps++;
00939     }
00940 
00941     // if (m_statushit == CRACK) m_cracks++;
00942     if (m_status == EMPTY) m_gaps++;
00943 
00944     if (m_status == FOUND && numDataPoints() > 0) {
00945         m_lstGaps =0;
00946 
00947         m_lstLayer = lastKPlane().getIDPlane();
00948         int indexHit = lastKPlane().getIDHit();
00949         int type = GFtutor::_DATA->clusterNoise(m_axis, indexHit);
00950         if(type > 0) {
00951             m_noisyHits++;
00952             if (kplanelist.size() < 3 ) m_noisyHits++;
00953         }
00954     }
00955 
00956     if (m_status == FOUND && numDataPoints() == 0) CONTROL_error++;
00957 
00958 }
00959 
00960 //#########################################################################
00961 void GFtrack::fit()
00962 //#########################################################################
00963 {
00964 
00965     GFdata::ini();
00966 
00967     if (kplanelist.size()< 3) {
00968         KalTrack::clear();
00969         return;
00970     }
00971 
00972     //----------------------------
00973     // Voila monsieur Kalman
00974     //----------------------------
00975     KalTrack::doFit();
00976     //----------------------------
00977 
00978     loadGFdata();
00979 }
00980 //#########################################################################
00981 void GFtrack::loadGFdata()
00982 //#########################################################################
00983 {
00984 
00985     // m_qbest = doQbest();
00986 
00987     // Output Data
00988     double x0,y0,z0;
00989     double slopex,slopey;
00990     x0=y0=z0=0.;
00991     slopex=slopey=0.;
00992     z0=kplanelist[0].getZPlane();
00993     if (m_axis == SiData::X) {
00994         x0=position(0.);
00995         slopex=slope();
00996     } else {
00997         y0=position(0.);
00998         slopey=slope();
00999     }
01000     m_vertex = Point(x0,y0,z0);
01001     double factor = -1.;
01002     m_direction = Vector(factor*slopex,factor*slopey,factor).unit();
01003     m_RCenergy =  KalEnergy();
01004     m_quality = computeQuality();
01005     m_firstLayer = kplanelist[0].getIDPlane();
01006     m_nhits = numDataPoints();
01007     m_itower = kplanelist[0].getIDTower();
01008 
01009 }
01010 
01011 //#########################################################################
01012 double GFtrack::computeQuality() const
01013 //#########################################################################
01014 {
01015     double quality = 0.;
01016     quality = 27./(9.+chiSquare()) + 15./(5.+chiSquareSegment()) +
01017         2.*(numDataPoints()-3.) - numGaps() - 5.*numFirstGaps();
01018     //          + (1./(0.1+abs(kink(0))));
01019     return quality;
01020 }
01021 //#########################################################################
01022 bool GFtrack::end() const
01023 //#########################################################################
01024 {
01025     bool end = !alive();
01026     if (m_lstGaps >= CONTROL_maxConsecutiveGaps) end = true;
01027     return end;
01028 }               
01029 
01030 //#########################################################################
01031 void GFtrack::kill()
01032 //#########################################################################
01033 {
01034     if (m_status == EMPTY) {
01035         m_gaps+=(-m_lstGaps);
01036         m_lstGaps = 0;
01037     }
01038     m_status = EMPTY;
01039 
01040     m_alive = false;
01041 
01042 }
01043 //########################################################
01044 void GFtrack::setAlive()
01045 //########################################################
01046 {
01047     m_alive = true;
01048 }
01049 //########################################################
01050 void GFtrack::setIniEnergy(double ene)
01051 //########################################################
01052 {
01053     m_iniEnergy = ene;
01054     KalTrack::setIniEnergy(ene);
01055 }
01056 
01057 //################################################
01058 KalPlane GFtrack::firstKPlane() const
01059 //################################################
01060 {
01061     if (kplanelist.size() == 0) {
01062         std::cout << "ERROR GFtrack::thisKPlane " << "\n";
01063         return originalKPlane();
01064     }
01065     return kplanelist.front();
01066 }
01067 
01068 //################################################
01069 KalPlane GFtrack::lastKPlane() const
01070 //################################################
01071 {
01072     if (kplanelist.size() == 0) {
01073         return originalKPlane();
01074     }
01075     return kplanelist.back();
01076 }
01077 //################################################
01078 KalPlane GFtrack::previousKPlane() const
01079 //################################################
01080 {
01081     if (kplanelist.size() <= 1) {
01082         //              std::cout << "ERROR GFtrack::previousKPlane " << "\n";
01083         return originalKPlane();
01084     }
01085     int iprevious = kplanelist.size()-2;
01086     if (iprevious == -1) return originalKPlane();
01087     return kplanelist[iprevious];
01088 }
01089 //################################################
01090 KalPlane GFtrack::originalKPlane() const
01091 //################################################
01092 {
01093 
01094     Ray testRay(m_inVertex,m_inDirection);
01095 
01096     double x0=testRay.position(0.).x();
01097     double y0=testRay.position(0.).y();
01098     double z0=testRay.position(0.).z();
01099 
01100     double dirX=testRay.direction().x();
01101     double dirY=testRay.direction().y();
01102     double dirZ=testRay.direction().z();
01103 
01104     double x = 0.;
01105     double x_orth = 0.;
01106     m_axis == SiData::X? x = x0: x=y0;
01107     m_axis == SiData::X? x_orth = y0: x_orth = x0;
01108 
01109     double slope=0.;
01110     double slope_orth=0.;
01111     if (dirZ != 0.) {
01112         slope = (m_axis == SiData::X? dirX/dirZ: dirY/dirZ);
01113         slope_orth = (m_axis == SiData::X? dirY/dirZ: dirX/dirZ);
01114     }
01115 
01116     KalPar porth(x_orth, slope_orth);
01117 
01118     KalPar pfit(x,slope);
01119 
01120     //unused:    double sigma = SiTracker::siStripPitch()/sqrt(12.);
01121     /*KalMatrix Q=KalPlane::Q(m_iniEnergy,slope,slope_orth,
01122     KalPlane::radLen(m_iniLayer)); */
01123     double sigma2Slope = CONTROL_iniErrorSlope * CONTROL_iniErrorSlope;
01124     double sigma2Position = CONTROL_iniErrorPosition * CONTROL_iniErrorPosition;
01125     KalMatrix covfit(sigma2Position,sigma2Slope,0.);
01126 
01127     KalHit hitfit(KalHit::FIT, pfit, covfit);
01128 
01129     KalPlane kp(0,m_iniLayer-1,m_iniEnergy, z0, porth, hitfit);
01130 
01131     return kp;
01132 }
01133 
01134 //#######################################################
01135 void GFtrack::removeStep(int iplane)
01136 //#######################################################
01137 {
01138     if (iplane == -1) iplane = numDataPoints()-1;
01139     if (iplane == -1) CONTROL_error++;
01140 
01141     if (iplane == numDataPoints()-1) {
01142         kplanelist.pop_back();
01143     } else {
01144         // WORK : remove the plane k
01145     }
01146 
01147     setStatus(EMPTY);
01148 }
01149 
01150 //#########################################################
01151 double GFtrack::doQbest()
01152 //#########################################################
01153 {
01154     double qbest = -1e6;
01155     if (numDataPoints()<3) return qbest;
01156     KalTrack::doFit();
01157     loadGFdata();
01158     // qbest=-1.*chiSquareSegment(m_sigmaCut*m_sigmaCut);
01159     m_qbest = Q();
01160     GFdata::ini();
01161     return m_qbest;
01162 }
01163 
01164 //#########################################################################
01165 void GFtrack::associateOrthStep(const GFtrack* _GFtrk, KalHit::TYPE typ)
01166 //#########################################################################
01167 {
01168     // locate the FIT hit of the previous plane in this step!
01169     if (numDataPoints() >=1 && _GFtrk->numDataPoints() >=1 && m_status == FOUND) {
01170 
01171         KalPar XPar = _GFtrk->lastKPlane().getHit(typ).getPar();
01172 
01173         double dz = kplanelist.back().getZPlane()-
01174             _GFtrk->kplanelist.back().getZPlane();
01175         double x0 = XPar.getPosition()+dz*XPar.getSlope();
01176         KalPar OrthPar(x0,XPar.getSlope());
01177 
01178         kplanelist.back().setOrthPar(OrthPar); // it changes the plane contents
01179     }
01180 }
01181 
01182 //#########################################################################
01183 void GFtrack::associateOrthGFtrack(const GFtrack* _GFtrk, bool fix, KalHit::TYPE typ)
01184 //#########################################################################
01185 {
01186     if (numDataPoints() == 0 || _GFtrk->numDataPoints() == 0) return;
01187     // Associate another track
01188     double XLIMIT = -91.;
01189     double SLOPENULL = 0.;
01190     int mplanes = _GFtrk->numDataPoints();
01191     for ( int iplane = 0; iplane < numDataPoints(); iplane++) {
01192         int idplane = kplanelist[iplane].getIDPlane();
01193         int jplane = -1;
01194         int jdplane = 0;
01195         do {
01196             jplane++;
01197             jdplane = _GFtrk->kplanelist[jplane].getIDPlane();
01198         } while (jdplane <= idplane && jplane+1 < mplanes-1);
01199         if (jdplane > idplane && jplane > 0) jplane--;
01200         
01201         
01202         KalPar Ypar = _GFtrk->kplanelist[jplane].getHit(typ).getPar();
01203         double dz = kplanelist[iplane].getZPlane()
01204             - _GFtrk->kplanelist[jplane].getZPlane();
01205         double x0 = Ypar.getPosition()+dz*Ypar.getSlope();
01206         double slope0 = Ypar.getSlope();
01207         if (!fix) {
01208             x0 = XLIMIT;
01209             slope0 = SLOPENULL;
01210         }
01211         KalPar OrthPar(x0,slope0);
01212 
01213         kplanelist[iplane].setOrthPar(OrthPar);
01214     }
01215 }
01216 //-----------------------------------------------------
01217 //
01218 //   GFparticle
01219 //
01220 //-----------------------------------------------------
01221 //######################################################
01222 GFparticle::GFparticle(double sigmaCut,
01223                        double energy,
01224                        int ist,
01225                        const Ray& testRay,
01226                        bool run) : GFbase(sigmaCut,energy,ist,testRay)
01227                        //#######################################################
01228 {
01229     ini();
01230     if (run == true) {
01231         doit();
01232         fit();
01233         if (empty()) clear();
01234     }
01235 }
01236 
01237 //##########################################
01238 void GFparticle::flagAllHits()
01239 //##########################################
01240 {
01241     _mXGFtrack->flagAllHits();
01242     _mYGFtrack->flagAllHits();
01243 }
01244 
01245 //##########################################
01246 void GFparticle::unFlagAllHits()
01247 //##########################################
01248 {
01249     _mXGFtrack->unFlagAllHits();
01250     _mYGFtrack->unFlagAllHits();
01251 }
01252 //########################################################
01253 bool GFparticle::empty() const
01254 //########################################################
01255 {
01256     bool empty = GFdata::empty();
01257     if (empty) return empty;
01258     empty = _mXGFtrack->empty() || _mYGFtrack->empty();
01259     return empty;
01260 }
01261 
01262 //########################################################
01263 bool GFparticle::accept() const
01264 //########################################################
01265 {
01266     bool valid = false;
01267     if (empty()) return valid;
01268 
01269     if (Q() < CUT_minQuality) return valid;
01270 
01271     if (GFtutor::CUT_veto) {
01272         int idhitX = -1;
01273         double sigmaX = -1.;
01274         int idhitY = -1;
01275         double sigmaY = -1.;
01276         if (_mXGFtrack->veto(idhitX,sigmaX) &&
01277             _mYGFtrack->veto(idhitY,sigmaY)) return valid;
01278     }
01279 
01280     valid = true;
01281     return valid;
01282 }
01283 //##########################################
01284 void GFparticle::clear()
01285 //##########################################
01286 {
01287     _mXGFtrack->clear();
01288     _mYGFtrack->clear();
01289 
01290     m_gaps   = 0;
01291     m_istGaps= 0;
01292     m_noisyHits = 0;
01293     m_istNoisyHits   = 0;
01294     m_lstLayer = 0;
01295 
01296     m_status = EMPTY;
01297     m_associate = true;
01298     m_conflictPattern = false;
01299 
01300     m_qbest = -1e6;
01301     GFdata::ini();
01302     setAlive();
01303 
01304 }
01305 //########################################################
01306 void GFparticle::writeOut(int level) const
01307 //########################################################
01308 {
01309     // kludge to avoid egcs warning messages
01310     int status = m_status;
01311     std::cout << " --- GFparticle::writeOut --- " << "\n";
01312     std::cout << " Qbest          = " << Qbest() << "\n";
01313     std::cout << " last  Layer    = " << lastLayer() << "\n";
01314     std::cout << " num Gaps       = " << numGaps() << "\n";
01315     std::cout << " num First Gaps = " << numFirstGaps() << "\n";
01316     std::cout << " num Noise      = " << numNoise() << "\n";
01317     std::cout << " num First Noise= " << numFirstNoise() << "\n";
01318     std::cout << " last Status    = " << status << "\n";
01319 
01320     GFdata::writeOut(level);
01321 
01322     std::cout << " -->  X - GFtrack : " << "\n";
01323     _mXGFtrack->writeOut(level);
01324     std::cout << " -->  Y - GFtrack : " << "\n";
01325     _mYGFtrack->writeOut(level);
01326 }
01327 
01328 //--------------------------------------------------------
01329 //  GFparticle - Private
01330 //--------------------------------------------------------
01331 
01332 //########################################################
01333 void GFparticle::ini()
01334 //########################################################
01335 {
01336     m_gaps   = 0;
01337     m_istGaps= 0;
01338     m_noisyHits = 0;
01339     m_istNoisyHits   = 0;
01340     m_lstLayer = 0;
01341 
01342     m_status = EMPTY;
01343     m_associate = true;
01344     m_conflictPattern = false;
01345 
01346     Ray testRay(m_inVertex,m_inDirection);      
01347     _mXGFtrack = new GFtrack(SiData::X, m_sigmaCut,
01348         m_iniEnergy, m_iniLayer, testRay, false);
01349 
01350     _mYGFtrack = new GFtrack(SiData::Y, m_sigmaCut,
01351         m_iniEnergy, m_iniLayer, testRay, false);
01352 
01353     m_qbest = -1e6;
01354     GFdata::ini();
01355     setAlive();
01356 
01357 }
01358 
01359 //#######################################################
01360 void GFparticle::step(int kplane)
01361 //#######################################################
01362 {
01363 
01364     if (!alive()) return;
01365 
01366     _mXGFtrack->step(kplane);
01367     _mYGFtrack->step(kplane);
01368 
01369     if (m_associate) {          
01370         associateStatus();
01371         associateStep();
01372     }
01373 }
01374 //#########################################################################
01375 void GFparticle::anastep(int kplane)
01376 //#########################################################################
01377 {
01378     if (!alive()) return;
01379 
01380     _mXGFtrack->anastep(kplane);
01381     _mYGFtrack->anastep(kplane);
01382 
01383     contability(kplane);
01384 
01385     if (m_associate) {
01386         associateAnaStep();
01387     }
01388 
01389     if (end()) {
01390         kill();
01391     }
01392 }
01393 
01394 //#########################################################################
01395 void GFparticle::fit()
01396 //#########################################################################
01397 {
01398     GFdata::ini();
01399 
01400     _mXGFtrack->fit();
01401     _mYGFtrack->fit();
01402 
01403     if (!_mXGFtrack->empty() && !_mYGFtrack->empty())
01404         loadGFdata();
01405 
01406     if (m_associate &&
01407         (!_mXGFtrack->empty() && !_mYGFtrack->empty()))
01408         associateFit();
01409 }
01410 
01411 //#########################################################################
01412 bool GFparticle::end() const
01413 //#########################################################################
01414 {
01415     bool end = !alive();
01416     // This comparatios is between different status and makes no sense
01417     // if (m_status == DONE) return end = true;
01418     if (!_mXGFtrack->alive() || !_mYGFtrack->alive()) end = true;
01419     return end;
01420 }
01421 
01422 //#########################################################################
01423 void GFparticle::kill()
01424 //#########################################################################
01425 {
01426     m_alive = false;
01427     _mXGFtrack->kill();
01428     _mYGFtrack->kill();
01429 
01430 }
01431 //#########################################################################
01432 void GFparticle::setAlive()
01433 //#########################################################################
01434 {
01435     m_alive = true;
01436     _mXGFtrack->setAlive();
01437     _mYGFtrack->setAlive();
01438 
01439 }
01440 
01441 
01442 //#########################################################################
01443 void GFparticle::loadGFdata()
01444 //#########################################################################
01445 {
01446     int ixlayer = _mXGFtrack->firstLayer();
01447     int iylayer = _mYGFtrack->firstLayer();
01448 
01449     m_firstLayer = (ixlayer <= iylayer? ixlayer : iylayer);
01450 
01451     int ixtower = _mXGFtrack->tower();
01452     //    int iytower = _mYGFtrack->tower();
01453     m_itower = ixtower;
01454 
01455     m_nhits = _mXGFtrack->nhits() + _mYGFtrack->nhits();
01456 
01457     m_quality = _mXGFtrack->Q() + _mYGFtrack->Q();
01458 
01459     m_RCenergy = 0.5*(_mXGFtrack->RCenergy() + _mYGFtrack->RCenergy());
01460 
01461     Ray XRay = _mXGFtrack->ray();
01462     Ray YRay = _mYGFtrack->ray();
01463     m_vertex=GFbase::doVertex(XRay, YRay);
01464 
01465     m_direction = GFbase::doDirection(_mXGFtrack->direction(),_mYGFtrack->direction());
01466 }
01467 //#########################################################################
01468 void GFparticle::contability(int kplane)
01469 //#########################################################################
01470 {
01471     if (!alive()) return;
01472 }
01473 
01474 //#########################################################################
01475 void GFparticle::associateStep()
01476 //#########################################################################
01477 {
01478     bool ok = false;
01479     bool done = false;
01480 
01481     ok = GFparticle::sameTower(_mXGFtrack,_mYGFtrack);
01482     if (!ok) done = GFparticle::removeWorseStep(_mXGFtrack,_mYGFtrack);
01483 
01484     if (!ok && !done) m_conflictPattern = true;
01485 
01486 }
01487 
01488 //#########################################################################
01489 void GFparticle::associateStatus()
01490 //#########################################################################
01491 {
01492     if (_mXGFtrack->status() == FOUND || _mYGFtrack->status() == FOUND)
01493         m_status = FOUND;
01494     else m_status = EMPTY;
01495 }
01496 
01497 //#########################################################################
01498 void GFparticle::associateAnaStep()
01499 //#########################################################################
01500 {
01501     _mXGFtrack->associateOrthStep(_mYGFtrack);  
01502     _mYGFtrack->associateOrthStep(_mXGFtrack);
01503 }
01504 
01505 //#########################################################################
01506 void GFparticle::associateFit()
01507 //#########################################################################
01508 {
01509     _mXGFtrack->associateOrthGFtrack(_mYGFtrack, m_associate, KalHit::SMOOTH);  
01510     _mYGFtrack->associateOrthGFtrack(_mXGFtrack, m_associate, KalHit::SMOOTH);
01511 }
01512 
01513 //#########################################################################
01514 bool GFparticle::sameTower(const GFtrack* _GFtrack1, const GFtrack* _GFtrack2)
01515 //#########################################################################
01516 {
01517     bool sametower = true;
01518     int itower1 = 0;
01519 
01520     if (_GFtrack1->status() == FOUND && _GFtrack1->numDataPoints() > 0) {
01521         itower1 = _GFtrack1->lastKPlane().getIDTower();
01522     }
01523     int itower2 = 0;
01524 
01525     if (_GFtrack2->status() == FOUND && _GFtrack2->numDataPoints() > 0) {
01526         itower2 = _GFtrack2->lastKPlane().getIDTower();
01527     }
01528 
01529     if (itower1 >= 10 && itower2 >= 10 && itower1 != itower2) sametower = false;
01530     return sametower;
01531 }
01532 //#########################################################################
01533 bool GFparticle::removeWorseStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2)
01534 //#########################################################################
01535 {
01536     bool done = false;
01537 
01538     if (_GFtrack1->status() != FOUND || _GFtrack2->status() != FOUND) CONTROL_error++;
01539     if (_GFtrack1->numDataPoints() == 0 || _GFtrack2->numDataPoints() == 0) CONTROL_error++;
01540 
01541     done = true;
01542 
01543     int kplane1 = _GFtrack1->previousKPlane().getIDPlane();
01544     int kplane2 = _GFtrack2->previousKPlane().getIDPlane();
01545 
01546     // the closest track first
01547     if (kplane1 != kplane2) {
01548         if (kplane1 > kplane2) _GFtrack2->removeStep();
01549         else  _GFtrack1->removeStep();
01550         return done;
01551     }
01552 
01553     int nhits1 = _GFtrack1->numDataPoints();
01554     int nhits2 = _GFtrack2->numDataPoints();
01555 
01556     double chi1 = _GFtrack1->lastKPlane().getDeltaChiSq(KalHit::FIT);
01557     double chi2 = _GFtrack2->lastKPlane().getDeltaChiSq(KalHit::FIT);
01558     if (nhits1 < 2 || nhits2 < 2) {
01559         if (_GFtrack1->_mGFsegment->accept() ||
01560             _GFtrack2->_mGFsegment->accept()  ) {
01561             if (_GFtrack1->_mGFsegment->accept()) chi1 = _GFtrack1->_mGFsegment->chiGFSq();
01562             else chi1 = 1e6;
01563             if (_GFtrack2->_mGFsegment->accept()) chi2 = _GFtrack2->_mGFsegment->chiGFSq();
01564             else chi2 = 1e6;
01565         }
01566     }
01567 
01568     if (chi2 > chi1) _GFtrack2->removeStep();
01569     else if (chi2 < chi1 ) _GFtrack1->removeStep();
01570     else done = false;
01571 
01572     return done;
01573 }
01574 //-------------------------------------------------------------------------
01575 //
01576 //     GFpair
01577 //
01578 //-------------------------------------------------------------------------
01579 
01580 //#########################################################################
01581 GFpair::GFpair(double xene, enum SiData::Axis axis,
01582                double sigmaCut,
01583                double energy,
01584                int ist,
01585                const Ray& testRay,
01586                bool run) : GFbase(sigmaCut,energy,ist,testRay),
01587                m_xEne(xene),
01588                m_axis(axis)
01589 //##########################################################################
01590 {
01591     ini();
01592     if (run == true) {
01593         doit();
01594         fit();
01595         if (empty()) clear();
01596     }
01597 }
01598 //#########################################################################
01599 void GFpair::flagAllHits()
01600 //#########################################################################
01601 {
01602     _mGFbest->flagAllHits();
01603     _mGFpair->flagAllHits();
01604 }
01605 //#########################################################################
01606 void GFpair::unFlagAllHits()
01607 //#########################################################################
01608 {
01609     _mGFbest->unFlagAllHits();
01610     _mGFpair->unFlagAllHits();
01611 }
01612 //#########################################################################
01613 bool GFpair::empty() const
01614 //#########################################################################
01615 {
01616     bool empty = GFdata::empty();
01617     if (empty) return empty;
01618 
01619     empty = empty || _mGFbest->empty();
01620     return empty;
01621 }
01622 
01623 //#########################################################################
01624 bool GFpair::accept() const
01625 //#########################################################################
01626 {
01627     bool ok = false;
01628     if (_mGFbest->empty()) return ok;
01629 
01630     // if there is a pair please do not impose the veto!
01631     int indexhit;
01632     double sigma;
01633     ok = true;
01634     if (GFtutor::CUT_veto) ok = !_mGFbest->veto(indexhit,sigma);
01635 
01636     return ok;
01637 }
01638 
01639 //#########################################################################
01640 void GFpair::clear()
01641 //#########################################################################
01642 {
01643     m_status = TOGETHER;
01644     m_decideBest = false;
01645 
01646     m_weightBest = 0.;
01647     m_errorSlope = 0.;
01648 
01649     m_together = 0;
01650     m_split = 0;
01651     m_one = 0;
01652     m_shared = 0;
01653     m_empty = 0;
01654 
01655     _mGFbest->clear();
01656     _mGFpair->clear();  
01657     if (_mGFalive) _mGFalive->clear();
01658 
01659     GFdata::ini();
01660     setAlive();
01661 }
01662 
01663 //########################################################
01664 void GFpair::writeOut(int level) const
01665 //########################################################
01666 {
01667     // kludge to avoid warnings from egcs
01668     int axis   = m_axis;
01669     int status = m_status;
01670     std::cout << " --- GFpair::writeOut --- " << "\n";
01671     std::cout << " axis:           " << axis << "\n";
01672     std::cout << " Energy Split    " << m_xEne << "\n";
01673     std::cout << " Weight          " << m_weightBest << "\n";
01674     std::cout << " planes together " << numTogether() << "\n";
01675     std::cout << " planes split    " << numSplit() << "\n";
01676     std::cout << " planes one      " << numOne() << "\n";
01677     std::cout << " shared hits     " << numSharedHits() << "\n";
01678     std::cout << " empty planes    " << numEmpty() << "\n";
01679     std::cout << " last Status     " << status << "\n";
01680 
01681     GFdata::writeOut(level);
01682     std::cout << " --> best Track : " << "\n";
01683     _mGFbest->writeOut(level);
01684     std::cout << " --> pair Track : " << "\n";
01685     _mGFpair->writeOut(level);
01686 }
01687 //-------------------------------------------------------------------------
01688 //    GFpair - Private
01689 //-------------------------------------------------------------------------
01690 //#########################################################################
01691 void GFpair::ini()
01692 //#########################################################################
01693 {
01694     // status
01695     m_status = TOGETHER;
01696     m_decideBest = false;
01697 
01698     // contability
01699     m_together = 0;
01700     m_split = 0;
01701     m_one = 0;
01702     m_shared = 0;
01703     m_empty = 0;
01704 
01705     // pair addition results
01706     m_weightBest = 0.;
01707     m_errorSlope = 0.;
01708 
01709     // internal variables
01710     _mGFbest = 0;
01711     _mGFpair = 0;
01712     _mGFalive = 0;
01713 
01714     Ray testRay(m_inVertex,m_inDirection);      
01715 
01716     _mGFbest = new GFtrack(m_axis, CONTROL_maxSigmaCut,
01717         CONTROL_xEne*m_iniEnergy, m_iniLayer, testRay, false);
01718 
01719     _mGFpair = new GFtrack(m_axis, CONTROL_maxSigmaCut,
01720         (1.-CONTROL_xEne)*m_iniEnergy, m_iniLayer, testRay, false);
01721 
01722     GFdata::ini();
01723     setAlive();
01724 }
01725 
01726 //#########################################################################
01727 void GFpair::step(int kplane)
01728 //#########################################################################
01729 {
01730     if (!alive()) return;
01731 
01732     newStatus(kplane);
01733 
01734     switch (status()) {
01735     case TOGETHER:
01736         stepTogether(kplane);
01737         break;
01738     case SPLIT:
01739         stepSplit(kplane);
01740         break;
01741     case ONE:
01742         _mGFalive->step(kplane);
01743         break;
01744     case DONE:
01745         break;
01746     }
01747 }
01748 
01749 //#########################################################################
01750 void GFpair::anastep(int kplane)
01751 //#########################################################################
01752 {       
01753     if (!m_alive) return;
01754 
01755     _mGFbest->anastep(kplane);
01756     _mGFpair->anastep(kplane);
01757 
01758     contability(kplane);
01759 
01760     if (end()) {
01761         kill();
01762     }
01763 }
01764 //#########################################################################
01765 void GFpair::fit()
01766 //#########################################################################
01767 {
01768     GFdata::ini();
01769     // Decide wich is the best track
01770     if (!m_decideBest) decideBest();
01771 
01772     if (m_decideBest) setIniEnergy();
01773     _mGFbest->fit();
01774     _mGFpair->fit();
01775 
01776     if (!_mGFbest->empty()) loadGFdata();
01777 }
01778 
01779 //#########################################################################
01780 bool GFpair::end() const
01781 //#########################################################################
01782 {
01783     bool end = !alive();
01784     if (m_status == DONE) return end = true;
01785     if (!_mGFbest->alive() && !_mGFpair->alive()) end = true;
01786     return end;
01787 }
01788 //#########################################################################
01789 void GFpair::kill()
01790 //#########################################################################
01791 {
01792     m_alive = false;
01793     _mGFbest->kill();
01794     _mGFpair->kill();
01795 }
01796 //#########################################################################
01797 void GFpair::setAlive()
01798 //#########################################################################
01799 {
01800     m_alive = true;
01801     _mGFbest->setAlive();
01802     _mGFpair->setAlive();
01803 }
01804 
01805 //#########################################################################
01806 void GFpair::contability(int kplane)
01807 //#########################################################################
01808 {
01809     if (m_status == TOGETHER) m_together++;
01810     else if (m_status == SPLIT) m_split++;
01811     else if (m_status == ONE) m_one++;
01812 
01813     if (m_status == SPLIT) {
01814         if (_mGFbest->status() == FOUND && _mGFpair->status() == FOUND) {
01815             if (_mGFbest->numDataPoints() == 0 || _mGFpair->numDataPoints() == 0) CONTROL_error++;
01816             else {
01817                 int idhit1 = _mGFbest->lastKPlane().getIDHit();
01818                 int idhit2 = _mGFpair->lastKPlane().getIDHit();
01819                 if (idhit1 == idhit2) m_shared++;
01820             }
01821         }
01822     }
01823 
01824     if (m_status == TOGETHER) {
01825         if (_mGFbest->status() != FOUND) m_empty++;
01826     } else if (m_status == ONE) {
01827         if (_mGFalive->status() != FOUND) m_empty++;
01828     } else if (m_status == SPLIT) {
01829         if (_mGFbest->status() != FOUND && _mGFpair->status() != FOUND) m_empty++;
01830     }
01831 
01832     if (m_status == ONE && m_split == 0 && m_together ==0) {
01833         m_one += m_together;
01834         m_together = 0;
01835     }
01836 
01837 }
01838 //#########################################################################
01839 void GFpair::loadGFdata()
01840 //#########################################################################
01841 {
01842     m_firstLayer = _mGFbest->firstLayer();
01843     m_itower = _mGFbest->tower();
01844     m_nhits = _mGFbest->nhits();
01845     if (!_mGFpair->empty()) m_nhits += _mGFpair->nhits()-m_shared;
01846 
01847     if (!_mGFpair->empty()) {   
01848         m_quality = _mGFbest->Q() + _mGFpair->Q()
01849             - (m_together/_mGFbest->numDataPoints())*(_mGFbest->Q());
01850     } else m_quality = _mGFbest->Q();
01851 
01852     m_RCenergy = doEnergy(_mGFbest, _mGFpair);
01853 
01854     if (_mGFpair->empty()) {
01855         m_direction = _mGFbest->direction();
01856         m_weightBest = 1.;
01857         m_errorSlope     = _mGFbest->errorSlopeAtVertex();
01858     } else {
01859         if (CONTROL_addTracksChi2) m_direction = doDirection(m_weightBest);
01860         else  m_direction = doDirection(_mGFbest, _mGFpair, m_weightBest, m_errorSlope);
01861         // m_direction = doDirectionXene(m_RCenergy, m_weightBest);
01862     }
01863 
01864 
01865     m_vertex = _mGFbest->vertex();
01866     if (!_mGFpair->empty()) {   
01867         Ray bestRay = Ray(_mGFbest->vertex(),_mGFbest->direction());
01868         Ray pairRay = Ray(_mGFpair->vertex(),_mGFpair->direction());
01869         m_vertex=GFbase::doVertex(bestRay, pairRay);
01870     }
01871 }
01872 //########################################################################
01873 void GFpair::newStatus(int kplane)
01874 //#########################################################################
01875 {
01876     StatusPair newStatus = status();
01877 
01878     switch(status()){
01879     case TOGETHER:
01880         if (!_mGFbest->alive()) newStatus = DONE;
01881         else if (forceSplit(kplane)) newStatus = ONE;
01882         break;
01883     case SPLIT:
01884         // You can be split but with only one hit, the together step should be used.
01885         if (_mGFbest->numDataPoints() <= 0) newStatus = TOGETHER;
01886         if (!_mGFbest->alive() && !_mGFpair->alive()) newStatus = DONE;
01887         else if (!_mGFbest->alive() || !_mGFpair->alive()) newStatus = ONE;
01888         break;
01889     case ONE:
01890         if (!_mGFalive->alive()) newStatus = DONE;
01891         break;
01892     case DONE:
01893         break;
01894     }
01895 
01896     setStatus(newStatus);
01897 }
01898 //########################################################################
01899 bool GFpair::forceSplit(int kplane) const
01900 //########################################################################
01901 {
01902     // WORK : define the best spliting point according with the energy
01903     bool split = false;
01904     if ( kplane > CONTROL_minSegmentHits + m_iniLayer -1) split = true;
01905     return split;
01906 }
01907 //########################################################################
01908 void GFpair::setStatus(StatusPair newStatus)
01909 //########################################################################
01910 {
01911     if (newStatus == status()) return;
01912 
01913     if (newStatus == ONE) {
01914         if (m_status == TOGETHER) {
01915             _mGFpair->clear();
01916             _mGFpair->kill();
01917         }
01918         _mGFalive = (_mGFbest->alive()? _mGFbest : _mGFpair);
01919     } else if (newStatus == DONE) {
01920         kill();
01921     }
01922 
01923     if (m_status == DONE || newStatus  == DONE) m_status = DONE;
01924     else if (m_status == ONE || newStatus == ONE) m_status = ONE;
01925     else if (m_status == SPLIT || newStatus == SPLIT) {
01926         m_status = SPLIT;
01927         if (_mGFbest->numDataPoints() == 0) m_status = TOGETHER;
01928     }
01929     else m_status = TOGETHER;
01930 
01931     // error control
01932     if (m_status == ONE) {
01933         if (_mGFbest->alive() && _mGFpair->alive()) CONTROL_error++;
01934     } else if (m_status != DONE) {
01935         if (!_mGFbest->alive() || !_mGFpair->alive()) CONTROL_error++;
01936     } else {
01937         if (_mGFbest->alive() || _mGFpair->alive()) CONTROL_error++;
01938     }
01939 }
01940 //#########################################################################
01941 void GFpair::stepTogether(int klayer)
01942 //#########################################################################
01943 {
01944     _mGFbest->step(klayer);
01945     if (!_mGFbest->_mGFsegment->accept()) return;
01946 
01947     _mGFbest->_mGFsegment->flagUsedHits(klayer);
01948     _mGFpair->_mGFsegment->best(klayer);
01949     _mGFbest->_mGFsegment->unFlagAllHits();
01950 
01951     bool split = _mGFpair->_mGFsegment->accept();
01952     if (split && _mGFpair->numDataPoints() == 0 ) {
01953         if (_mGFpair->_mGFsegment->getKPlane().getSigma(KalHit::SMOOTH) > sigmaCut()) {
01954             split= false;
01955         }
01956     }
01957 
01958     if (split) {
01959         setStatus(SPLIT);
01960         _mGFpair->kplanelist.push_back(_mGFpair->_mGFsegment->getKPlane());
01961     } else {
01962         _mGFpair->step(klayer);
01963     }
01964     _mGFpair->setStatus(FOUND);
01965 
01966 }
01967 //#########################################################################
01968 void GFpair::stepSplit(int klayer)
01969 //#########################################################################
01970 {
01971     _mGFbest->step(klayer);
01972     _mGFpair->step(klayer);
01973 
01974     if (_mGFbest->status() != FOUND || _mGFpair->status() != FOUND) return;
01975 
01976     if (_mGFbest->numDataPoints()==0 || _mGFpair->numDataPoints() == 0) CONTROL_error++;
01977 
01978     int indexhit1 = _mGFbest->lastKPlane().getIDHit();
01979     int indexhit2 = _mGFpair->lastKPlane().getIDHit();
01980     if ( indexhit1 != indexhit2) return;
01981 
01982     selfishStepSplit(klayer);
01983 }
01984 
01985 //#########################################################################
01986 void GFpair::selfishStepSplit(int klayer)
01987 //#########################################################################
01988 {
01989 
01990     if (_mGFbest->numDataPoints() <2 || _mGFpair->numDataPoints()<2) CONTROL_error++;
01991 
01992     removeWorseStep(_mGFbest,_mGFpair);
01993     GFtrack* _GFwiner;
01994     GFtrack* _GFloser;
01995     if (_mGFpair->status() == EMPTY) {
01996         _GFwiner = _mGFbest;
01997         _GFloser = _mGFpair;
01998     } else {
01999         _GFwiner = _mGFpair;
02000         _GFloser = _mGFbest;
02001     }
02002 
02003     // flag the winer hit - remove the loser - next tray - unflag again
02004 
02005     int indexhit = _GFwiner->lastKPlane().getIDHit();
02006 
02007     GFtutor::_DATA->flagHit(m_axis, indexhit);
02008 
02009     _GFloser->step(klayer);
02010     GFtutor::_DATA->unflagHit(m_axis, indexhit);
02011 
02012     if (_GFloser->status() == EMPTY && allowedShareHit(_GFwiner)) {
02013         _GFloser->step(klayer);
02014         double sigma = _GFloser->lastKPlane().getSigma(KalHit::PRED);
02015         if (sigma > CUT_maxSigmasSharedHits*m_sigmaCut) {
02016             _GFloser->removeStep();
02017         }
02018     }
02019 
02020 }
02021 //#########################################################################
02022 void GFpair::removeWorseStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2)
02023 //#########################################################################
02024 {
02025     if (_GFtrack1->status() != FOUND || _GFtrack2->status() != FOUND) CONTROL_error++;
02026 
02027     //unused:   int indexhit1 = -1;
02028     //unused:   int indexhit2 = -1;
02029 
02030     int kplane1 = _GFtrack1->previousKPlane().getIDPlane();
02031     int kplane2 = _GFtrack2->previousKPlane().getIDPlane();
02032 
02033     // the closest track first
02034     if (kplane1 != kplane2) {
02035         if (kplane1 > kplane2) _GFtrack2->removeStep();
02036         else  _GFtrack1->removeStep();
02037         return;
02038     }
02039 
02040     int nhits1 = _GFtrack1->numDataPoints();
02041     int nhits2 = _GFtrack2->numDataPoints();
02042 
02043     double chi1 = _GFtrack1->lastKPlane().getDeltaChiSq(KalHit::FIT);
02044     double chi2 = _GFtrack2->lastKPlane().getDeltaChiSq(KalHit::FIT);
02045 
02046     if (nhits1 < 2 || nhits2 < 2) {
02047         if (_GFtrack1->_mGFsegment->accept() ||
02048             _GFtrack2->_mGFsegment->accept()  ) {
02049             if (_GFtrack1->_mGFsegment->accept()) chi1 = _GFtrack1->_mGFsegment->chiGFSq();
02050             else chi1 = 1e6;
02051             if (_GFtrack2->_mGFsegment->accept()) chi2 = _GFtrack2->_mGFsegment->chiGFSq();
02052             else chi2 = 1e6;
02053         }
02054     }
02055 
02056     if (chi2 >= chi1) {
02057         _GFtrack2->removeStep();
02058     } else {
02059         _GFtrack1->removeStep();
02060     }
02061 }
02062 //#########################################################################
02063 bool GFpair::allowedShareHit(const GFtrack* _GFtrack) const
02064 //#########################################################################
02065 {
02066     bool allow = false;
02067 
02068     if (_GFtrack->numDataPoints() == 0) CONTROL_error++;
02069 
02070     // the size of the cluster is big enough for this slope
02071     int idhit = _GFtrack->lastKPlane().getIDHit();
02072     double slope = _GFtrack->lastKPlane().getHit(KalHit::PRED).getPar().getSlope();
02073     int value = 2; // more strips than the expected
02074     if (GFtutor::okClusterSize(m_axis,idhit,slope) >= value) allow = true;
02075 
02076     return allow;
02077 }
02078 //#########################################################################
02079 Vector GFpair::doDirection(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2,
02080                            double& weight1, double& errorSlope)
02081 //#########################################################################
02082 {
02083     Vector dir(0.,0.,0.);
02084     weight1 = 0.;
02085     errorSlope = 0.;
02086     if (_GFtrk1->empty()) return dir;
02087 
02088     double xene = _GFtrk1->m_iniEnergy/(_GFtrk1->m_iniEnergy+_GFtrk2->m_iniEnergy);
02089     double slopeBest   = _GFtrk1->slope();
02090     double slopePair   = _GFtrk2->slope();
02091     double errorSQbest = _GFtrk1->errorSlopeAtVertex();
02092     errorSQbest *= errorSQbest;
02093     double errorSQpair = _GFtrk2->errorSlopeAtVertex();
02094     errorSQpair *= errorSQpair;
02095     double errorSQ     = errorSQbest*errorSQpair/((1.-xene)*errorSQbest+xene*errorSQpair);
02096 
02097     weight1     = xene*errorSQ/errorSQbest;
02098     double slope = errorSQ*(xene*(slopeBest/errorSQbest)+(1.-xene)*(slopePair/errorSQpair));
02099     errorSlope  = sqrt(errorSQ);
02100 
02101     double factor = -1;
02102     double slopex = 0;
02103     double slopey =0;
02104     if (_GFtrk1->getAxis() == SiData::X) slopex = slope;
02105     else slopey = slope;
02106 
02107     dir = Vector(factor*slopex,factor*slopey,factor).unit();
02108     return dir;
02109 }
02110 //########################################################
02111 Vector GFpair::doDirection(double& xFactor)
02112 //########################################################
02113 {
02114     Vector dir(0.,0.,0.);
02115     float wt1x = 1./(_mGFbest->chiSquare() + .01);
02116     float wt2x = 1./(3.*(_mGFpair->chiSquare() + .01));
02117     if (wt2x > wt1x) wt2x = wt1x;
02118     xFactor = wt1x/(wt1x + wt2x);
02119     Vector t1 = _mGFbest->direction();
02120     Vector t2 = _mGFpair->direction();
02121     //unused: float sinDLT12 = sqrt(std::max(0.0,(1. - sqr(t1*t2))));
02122     // if(sinDLT12 < .010/m_iniEnergy) { // only trust small opening angle paris
02123     double aveSlope = xFactor*_mGFbest->slope() + (1.-xFactor)*_mGFpair->slope();
02124     double slopeX = 0.;
02125     double slopeY = 0.;
02126     if (m_axis == SiData::X) slopeX = aveSlope;
02127     else slopeY = aveSlope;
02128     double factor = -1.;
02129     dir = Vector(factor*slopeX,factor*slopeY,factor).unit();
02130     // }  else dir = _mGFbest->direction();
02131     return dir;
02132 }
02133 //########################################################
02134 Vector GFpair::doDirectionXene(double xene, double& weight)
02135 //########################################################
02136 {
02137     Vector dir(0.,0.,0.);
02138     double x3 = xene*xene*xene;
02139     double ix3 = (1.-xene)*(1.-xene)*(1.-xene);
02140     float wt1x = x3/(x3+ix3);
02141     float wt2x = ix3/(x3+ix3);
02142     Vector t1 = _mGFbest->direction();
02143     Vector t2 = _mGFpair->direction();
02144     double aveSlope = wt1x*_mGFbest->slope() + wt2x*_mGFpair->slope();
02145     double slopeX = 0.;
02146     double slopeY = 0.;
02147     if (m_axis == SiData::X) slopeX = aveSlope;
02148     else slopeY = aveSlope;
02149     double factor = -1.;
02150     dir = Vector(factor*slopeX,factor*slopeY,factor).unit();
02151 
02152     weight = wt1x;
02153     return dir;
02154 }
02155 //#########################################################################
02156 double GFpair::doEnergy(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2)
02157 //#########################################################################
02158 {
02159     double ene1 = std::min((double)_GFtrk1->RCenergy(),m_iniEnergy);
02160     double ene2 = m_iniEnergy - ene1;
02161     if (!_GFtrk2->empty())
02162         ene2 = std::min((double)_GFtrk2->RCenergy(),m_iniEnergy);
02163     double x1 = ene1/m_iniEnergy;
02164     double x2 = 1.-ene2/m_iniEnergy;
02165     double x = 0.5*(x1+x2);
02166     return x;
02167 }
02168 
02169 //#########################################################################
02170 void GFpair::decideBest()
02171 //#########################################################################
02172 {
02173     // already decided!
02174     if (m_decideBest) return;
02175 
02176     double qbest = _mGFbest->doQbest();
02177     double qpair = _mGFpair->doQbest();
02178 
02179     if (qpair > qbest) swap();
02180 
02181     m_decideBest = true;
02182 }
02183 //#########################################################################
02184 void GFpair::swap()
02185 //#########################################################################
02186 {
02187 
02188     double eneBest = _mGFbest->m_iniEnergy;
02189     double enePair = _mGFpair->m_iniEnergy;
02190     double cutBest = _mGFbest->m_sigmaCut;
02191     double cutPair = _mGFbest->m_sigmaCut;
02192 
02193     _mGFalive = _mGFbest;
02194     _mGFbest  = _mGFpair;
02195     _mGFpair  = _mGFalive;
02196 
02197     _mGFbest->setIniEnergy(eneBest);
02198     _mGFpair->setIniEnergy(enePair);
02199     _mGFbest->m_sigmaCut = cutPair;
02200     _mGFpair->m_sigmaCut = cutBest;
02201 
02202     _mGFalive = (_mGFbest->m_alive? _mGFbest:_mGFpair);
02203 
02204 }
02205 
02206 //#########################################################################
02207 void GFpair::setIniEnergy()
02208 //#########################################################################
02209 {
02210     // set The real energies before the Fit!!
02211     _mGFbest->setIniEnergy(m_xEne*m_iniEnergy);
02212     _mGFpair->setIniEnergy((1.-m_xEne)*m_iniEnergy);
02213 }
02214 //-------------------------------------------------------------------------
02215 //
02216 //   Gamma
02217 //
02218 //-------------------------------------------------------------------------
02219 //#########################################################################
02220 GFgamma::GFgamma(double xene,
02221                  double sigmaCut,
02222                  double energy,
02223                  int ist,
02224                  const Ray& testRay) : GFbase(sigmaCut,energy,ist,testRay),
02225                  m_xEne(xene)
02226                  //##########################################################################
02227 {
02228 
02229     GFdata::ini();
02230 
02231     m_connect = GFtutor::CONTROL_connectGFpair;
02232     m_associate = false;
02233     m_patternSwap = false;
02234     construct();
02235     if (!m_conflictPattern) return;
02236 
02237     clear();
02238     delete _mXpair;
02239     delete _mYpair;
02240     _mXpair = 0;
02241     _mYpair = 0;
02242     m_associate = true;
02243     // m_patternSwap;
02244     construct();
02245 }
02246 //#########################################################################
02247 void GFgamma::flagAllHits()
02248 //#########################################################################
02249 {
02250     _mXpair->flagAllHits();
02251     _mYpair->flagAllHits();
02252 }
02253 //#########################################################################
02254 void GFgamma::unFlagAllHits()
02255 //#########################################################################
02256 {
02257     _mXpair->unFlagAllHits();
02258     _mYpair->unFlagAllHits();
02259 }
02260 //#########################################################################
02261 bool GFgamma::empty() const
02262 //##########################################################################
02263 {
02264     bool empty = GFdata::empty();
02265     empty = empty || _mXpair->empty() || _mYpair->empty();
02266     return empty;
02267 }
02268 
02269 //#########################################################################
02270 bool GFgamma::accept(const GFdata& pData1, const GFdata& pData2)
02271 //##########################################################################
02272 {
02273     bool ok = false;
02274 
02275     int ilayerXfirst = pData1.firstLayer();
02276     int ilayerYfirst = pData2.firstLayer();
02277     if (abs(ilayerXfirst-ilayerYfirst) > CONTROL_maxConsecutiveGaps) return ok;
02278 
02279     int iXtower = pData1.tower();
02280     int iYtower = pData2.tower();
02281     if (ilayerXfirst == ilayerYfirst) {
02282         if (iXtower == iYtower) ok = true;
02283         else ok = false;
02284     }
02285 
02286     return ok;
02287 }
02288 //#########################################################################
02289 void GFgamma::clear()
02290 //##########################################################################
02291 {
02292     m_status = TOGETHER;
02293     m_together = 0;
02294     m_split = 0;
02295     m_one = 0;
02296 
02297     _mXpair->clear();
02298     _mYpair->clear();
02299 
02300     GFdata::ini();
02301     setAlive();
02302 
02303 }
02304 
02305 
02306 //#########################################################################
02307 double GFgamma::Qbest()
02308 //##########################################################################
02309 {
02310     if (empty()) return m_quality;
02311     return _mXpair->_mGFbest->Qbest()+_mYpair->_mGFbest->Qbest();
02312 }
02313 //#########################################################################
02314 bool GFgamma::accept() const
02315 //##########################################################################
02316 {
02317     bool accept = false;
02318     if (empty()) return accept;
02319 
02320     accept = true;
02321     if (GFtutor::CUT_veto) accept = !veto();
02322 
02323     return accept;
02324 }
02325 
02326 //########################################################
02327 void GFgamma::writeOut(int level) const
02328 //########################################################
02329 {
02330     // kludge to avoid warning messages from egcs
02331     int status = m_status;
02332     std::cout << " --- GFgamma::writeOut --- " << "\n";
02333     std::cout << " planes together " << numTogether() <<"\n";
02334     std::cout << " planes split    " << numSplit() <<"\n";
02335     std::cout << " planes one      " << numOne() << "\n";
02336     std::cout << " last Status     " << status << "\n";
02337 
02338     GFdata::writeOut(level);
02339     std::cout << " --> Xpair : " <<"\n";
02340     _mXpair->writeOut(level);
02341     std::cout << " --> Ypair : " <<"\n";
02342     _mYpair->writeOut(level);
02343 
02344 }
02345 //#########################################################################
02346 bool GFgamma::veto() const
02347 //##########################################################################
02348 {
02349     bool veto = false;
02350     int ixhit;
02351     int iyhit;
02352 
02353     double sigma;
02354     // Veto for both views
02355     veto = _mXpair->_mGFbest->veto(ixhit,sigma) && _mYpair->_mGFbest->veto(iyhit,sigma);
02356     if (veto && SiClusters::tower(ixhit) != SiClusters::tower(iyhit)) veto = false;
02357 
02358     return veto;
02359 }
02360 
02361 //#########################################################################
02362 Point GFgamma::getFirstHit() const
02363 //##########################################################################
02364 {
02365     Point firstHit;
02366     double zx = _mXpair->_mGFbest->vertex().z();
02367     double zy = _mYpair->_mGFbest->vertex().z();
02368     double xx = _mXpair->_mGFbest->vertex().x();
02369     double yy = _mYpair->_mGFbest->vertex().y();
02370 
02371     if (zy > zx) xx = _mXpair->_mGFbest->position(zy-zx);
02372     else yy = _mYpair->_mGFbest->position(zx-zy);
02373 
02374     double zz = (zy > zx? zy : zx);
02375 
02376     firstHit = Point(xx,yy,zz);
02377     return firstHit;
02378 }
02379 
02380 //-------------------------------------------------------------------------
02381 //   Gamma - Private
02382 //-------------------------------------------------------------------------
02383 
02384 //#########################################################################
02385 void GFgamma::ini()
02386 //#########################################################################
02387 {
02388 
02389     _mXpair = 0;
02390     _mYpair = 0;
02391     _mXpair = new GFpair(m_xEne, SiData::X, m_sigmaCut,
02392         m_iniEnergy, m_iniLayer, Ray(m_inVertex,m_inDirection), false);
02393     _mYpair = new GFpair(m_xEne, SiData::Y, m_sigmaCut,
02394         m_iniEnergy, m_iniLayer, Ray(m_inVertex,m_inDirection), false);
02395 
02396     // controls
02397     m_fixTopology = false;
02398     m_decideBest = false;
02399     m_patternSwap = false;
02400     m_conflictPattern = false;
02401     m_swapDone = false;
02402 
02403     // status
02404     setDecideBest(m_decideBest);
02405     m_status = TOGETHER;
02406 
02407     //contability
02408     m_together = 0;
02409     m_split = 0;
02410     m_one = 0;
02411 
02412     setAlive();
02413     GFdata::ini();
02414 
02415 }
02416 
02417 //#########################################################################
02418 void GFgamma::step(int kplane)
02419 //#########################################################################
02420 {
02421     if (!m_alive) return;
02422 
02423     _mXpair->step(kplane);
02424     _mYpair->step(kplane);
02425 
02426     m_status = newStatus();
02427     if (m_connect) {
02428         connectStep();
02429         if (!m_fixTopology) topologyStep();
02430 
02431         if (m_associate) associateStep();
02432     }
02433 }
02434 
02435 //#########################################################################
02436 void GFgamma::anastep(int kplane)
02437 //#########################################################################
02438 {
02439     if (!m_alive) return;
02440 
02441     _mXpair->anastep(kplane);
02442     _mYpair->anastep(kplane);
02443 
02444     contability(kplane);
02445     if (m_associate) associateAnaStep();
02446     if (end()) {
02447         kill();
02448     }
02449 }
02450 
02451 //#########################################################################
02452 void GFgamma::fit()
02453 //##########################################################################
02454 {
02455     GFdata::ini();
02456 
02457     //  if (m_patternError) return;
02458     if (!m_decideBest) decideBest();
02459 
02460     _mXpair->fit();
02461     _mYpair->fit();
02462     if (_mYpair->empty()||_mXpair->empty()) return;
02463 
02464     loadGFdata();
02465     associateFit();
02466     //  if (m_associate) associateFit();
02467 }
02468 
02469 //#########################################################################
02470 bool GFgamma::end() const
02471 //#########################################################################
02472 {
02473     bool end = !m_alive;
02474     if (!_mXpair->m_alive || !_mYpair->m_alive) end = true;
02475     return end;
02476 }
02477 
02478 //#########################################################################
02479 void GFgamma::kill()
02480 //#########################################################################
02481 {
02482     m_alive = false;
02483     _mXpair->kill();
02484     _mYpair->kill();
02485 }
02486 
02487 //#########################################################################
02488 void GFgamma::setAlive()
02489 //#########################################################################
02490 {
02491     m_alive = true;
02492     _mXpair->setAlive();
02493     _mYpair->setAlive();
02494 }
02495 
02496 //#########################################################################
02497 void GFgamma::contability(int kplane)
02498 //#########################################################################
02499 {
02500     if (m_status == TOGETHER) m_together++;
02501     else if (m_status == SPLIT) m_split++;
02502     else if (m_status == ONE) m_one++;
02503 }
02504 
02505 //#########################################################################
02506 void GFgamma::loadGFdata()
02507 //##########################################################################
02508 {
02509     m_quality = _mXpair->Q() + _mYpair->Q();
02510 
02511     m_direction = GFdata::doDirection(_mXpair->direction(),_mYpair->direction());
02512 
02513     Ray XRay = Ray(_mXpair->vertex(),_mXpair->direction());
02514     Ray YRay = Ray(_mYpair->vertex(),_mYpair->direction());
02515     m_vertex=GFbase::doVertex( XRay, YRay);
02516 
02517     m_RCenergy = 0.5*(_mXpair->RCenergy()+_mYpair->RCenergy());
02518 
02519     int ixlayer = _mXpair->firstLayer();
02520     int iylayer = _mYpair->firstLayer();
02521     m_firstLayer = (ixlayer < iylayer? ixlayer:iylayer);
02522 
02523     m_nhits = _mXpair->nhits()+_mYpair->nhits();
02524 
02525     m_itower = ( ixlayer < iylayer? _mXpair->tower() : _mYpair->tower());
02526 }
02527 //#########################################################################
02528 void GFgamma::construct()
02529 //##########################################################################
02530 {
02531     ini();
02532 
02533     doit();
02534 
02535     fit();
02536 
02537     if (empty()) clear();
02538 }
02539 
02540 //#########################################################################
02541 GFbase::StatusPair GFgamma::newStatus()
02542 //#########################################################################
02543 {
02544     StatusPair Xstatus = _mXpair->status();
02545     StatusPair Ystatus = _mYpair->status();
02546     StatusPair status = Xstatus;
02547 
02548     if (Xstatus == Ystatus) status = Xstatus;
02549     else if (Xstatus == DONE || Ystatus == DONE) status = DONE;
02550     else if (Xstatus == SPLIT || Ystatus == SPLIT) status = SPLIT;
02551     else if (Xstatus == TOGETHER || Ystatus == TOGETHER) status = TOGETHER;
02552     return status;
02553 }
02554 //#########################################################################
02555 void GFgamma::connectStep()
02556 //#########################################################################
02557 {
02558     if (m_status == TOGETHER) {
02559 
02560         if (!GFparticle::sameTower(_mXpair->_mGFbest, _mYpair->_mGFbest)) {
02561 
02562             bool done = GFparticle::removeWorseStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
02563             if (!done) m_conflictPattern = true;
02564         }
02565     }
02566 
02567     if (!m_associate) {
02568         if (m_status == SPLIT) {
02569 
02570             bool bestbest = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFbest);
02571             bool pairpair = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFpair);
02572             if (!bestbest || !pairpair) {
02573 
02574                 bool bestpair = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFpair);
02575                 bool pairbest = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFbest);
02576                 if (!bestpair || !pairbest) m_conflictPattern = true;
02577             }
02578         }
02579         if (m_status == ONE) {
02580 
02581             if (!GFparticle::sameTower(_mXpair->_mGFalive,_mYpair->_mGFalive)) {
02582 
02583                 bool done = GFparticle::removeWorseStep(_mXpair->_mGFalive,_mYpair->_mGFalive);
02584                 if (!done) m_conflictPattern = true;
02585             }
02586         }
02587     }
02588 
02589 }
02590 
02591 //#########################################################################
02592 void GFgamma::topologyStep()
02593 //#########################################################################
02594 {
02595     if (m_status == SPLIT) {
02596         bool bestbest = crossingTowers(_mXpair->_mGFbest,_mYpair->_mGFbest,
02597             _mXpair->_mGFpair,_mYpair->_mGFpair);
02598         bool bestpair = crossingTowers(_mXpair->_mGFbest,_mYpair->_mGFpair,
02599             _mXpair->_mGFpair,_mYpair->_mGFbest);
02600         if (bestbest && bestpair) CONTROL_error++;
02601         else if (bestbest || bestpair) {
02602             m_fixTopology = true;
02603             if (!m_associate) {
02604                 m_associate = true;
02605                 m_patternSwap = bestpair;
02606             }
02607         }
02608     } else if (m_status == ONE) {
02609         if (_mXpair->_mGFalive->status() == FOUND &&
02610             _mYpair->_mGFalive->status() == FOUND) {
02611             m_fixTopology = true;
02612             if (!m_associate) {
02613                 m_associate = true;
02614                 if ((_mXpair->_mGFbest->alive() && _mYpair->_mGFbest->alive()) ||
02615                     (_mXpair->_mGFpair->alive() && _mYpair->_mGFpair->alive()) )
02616                     m_patternSwap = false;
02617                 else m_patternSwap = true;
02618             }
02619         }
02620     }
02621 }
02622 //#########################################################################
02623 void GFgamma::associateStep()
02624 //#########################################################################
02625 {
02626     associateStatus(m_status);
02627 
02628     if (m_status == SPLIT) {
02629         // TOGETHER already tested in connectStep;
02630         bool bestbest = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFbest);
02631         bool pairpair = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFpair);
02632         if (!bestbest || !pairpair) {
02633             bool done = false;
02634             // spetial case - first plane
02635 
02636             if (!bestbest) done = GFparticle::removeWorseStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
02637             if (!pairpair) done = GFparticle::removeWorseStep(_mXpair->_mGFpair,_mYpair->_mGFpair);
02638             if (!done) m_conflictPattern = true;
02639         }
02640     } else if (m_status == ONE) {
02641         if (!GFparticle::sameTower(_mXpair->_mGFalive, _mYpair->_mGFalive)) {
02642 
02643             bool done = GFparticle::removeWorseStep(_mXpair->_mGFalive,_mYpair->_mGFalive);
02644             if (!done) m_conflictPattern = true;
02645         }
02646     }
02647 }
02648 
02649 //#########################################################################
02650 void GFgamma::associateStatus(StatusPair status)
02651 //#########################################################################
02652 {       
02653     StatusPair Xstatus = _mXpair->status();
02654     StatusPair Ystatus = _mYpair->status();
02655     if (m_patternSwap && !m_swapDone && m_status != TOGETHER) {
02656         //      if (!_mXpair->_mGFpair->m_alive && !_mYpair->_mGFpair->m_alive) CONTROL_error++;
02657         if (Xstatus == SPLIT || _mXpair->numSplit()>0) _mXpair->swap();
02658         else if (Ystatus == SPLIT || _mYpair->numSplit()>0) _mYpair->swap();
02659         else CONTROL_error++;
02660         m_swapDone = true; // reset the pattern - now again best-best.
02661     }
02662 
02663     _mXpair->setStatus(status);
02664     _mYpair->setStatus(status);
02665 
02666     m_status = status;
02667 }
02668 //#########################################################################
02669 bool GFgamma::crossingTowers(const GFtrack* _GFtrkX1, const GFtrack* _GFtrkY1,
02670                              const GFtrack* _GFtrkX2, const GFtrack* _GFtrkY2)
02671                              //#########################################################################
02672 {
02673     // The two tracks are in different towers
02674     bool fix = false;
02675 
02676     int nX1hits = _GFtrkX1->numDataPoints();
02677     int nY1hits = _GFtrkY1->numDataPoints();
02678     if ( nX1hits < 1 || nY1hits < 1 ) return fix;
02679 
02680     int nX2hits = _GFtrkX2->numDataPoints();
02681     int nY2hits = _GFtrkY2->numDataPoints();
02682     if ( nX2hits < 1 || nY2hits < 1 ) return fix;
02683 
02684     if (_GFtrkX1->status() != FOUND || _GFtrkY1->status() != FOUND ||
02685         _GFtrkX2->status() != FOUND || _GFtrkY2->status() != FOUND ) return fix;
02686 
02687     int iX1tower = _GFtrkX1->kplanelist[nX1hits-1].getIDTower();
02688     int iY1tower = _GFtrkY1->kplanelist[nY1hits-1].getIDTower();
02689     int iX2tower = _GFtrkX2->kplanelist[nX2hits-1].getIDTower();
02690     int iY2tower = _GFtrkY2->kplanelist[nY2hits-1].getIDTower();
02691 
02692     if ((iX1tower == iY1tower) &&
02693         (iX2tower == iY2tower) &&
02694         (iX1tower != iX2tower)) fix = true;
02695 
02696     return fix;
02697 }
02698 
02699 
02700 //#########################################################################
02701 void GFgamma::associateAnaStep()
02702 //#########################################################################
02703 {
02704     associateAnaStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
02705     associateAnaStep(_mYpair->_mGFbest,_mXpair->_mGFbest);
02706     associateAnaStep(_mXpair->_mGFpair,_mYpair->_mGFpair);
02707     associateAnaStep(_mYpair->_mGFpair,_mXpair->_mGFpair);
02708 }
02709 
02710 //#########################################################################
02711 void GFgamma::associateAnaStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2)
02712 //#########################################################################
02713 {
02714     if (_GFtrack1->numDataPoints() == 0) return;
02715     if (_GFtrack1->status() != FOUND) return;
02716 
02717     if (!_GFtrack2->m_alive) {
02718         if (_GFtrack2->numDataPoints() < CONTROL_minSegmentHits) _GFtrack1->clear();
02719         _GFtrack1->kill();
02720     } else {
02721         _GFtrack1->associateOrthStep(_GFtrack2);
02722     }
02723 }
02724 //#########################################################################
02725 void GFgamma::associateFit()
02726 //#########################################################################
02727 {
02728     bool fix = m_fixTopology;
02729     KalHit::TYPE typ = KalHit::SMOOTH;
02730     _mXpair->_mGFbest->associateOrthGFtrack(_mYpair->_mGFbest, fix, typ);
02731     _mYpair->_mGFbest->associateOrthGFtrack(_mXpair->_mGFbest, fix, typ);
02732     _mXpair->_mGFpair->associateOrthGFtrack(_mYpair->_mGFpair, fix, typ);
02733     _mYpair->_mGFpair->associateOrthGFtrack(_mXpair->_mGFpair, fix, typ);
02734 }
02735 
02736 //#########################################################################
02737 void GFgamma::setDecideBest(bool decideBest)
02738 //##########################################################################
02739 {
02740     m_decideBest = decideBest;
02741     _mXpair->setDecideBest(m_decideBest);
02742     _mYpair->setDecideBest(m_decideBest);
02743 }
02744 //#########################################################################
02745 void GFgamma::decideBest()
02746 //##########################################################################
02747 {
02748 
02749     if (!m_associate) return;
02750 
02751     setDecideBest(true);
02752 
02753     double qbest = _mXpair->_mGFbest->doQbest();
02754     qbest += _mYpair->_mGFbest->doQbest();
02755 
02756     double qpair = _mXpair->_mGFpair->doQbest();
02757     qpair += _mYpair->_mGFpair->doQbest();
02758 
02759     if (qbest < qpair) {
02760         _mXpair->swap();
02761         _mYpair->swap();
02762     }
02763 
02764 }
02765 //--------------------------------------------------------
02766 //   MISCELANEOUS UTILITIES
02767 //--------------------------------------------------------
02768 //########################################################
02769 int GFtutor::okClusterSize(SiData::Axis axis, int indexhit, double slope)                       
02770 //########################################################
02771 {
02772     int icluster = 0;
02773 
02774     int size = (int) GFtutor::_DATA->clusterSize(axis,indexhit);
02775 
02776     double itemp =  si_strip_pitch;
02777 
02778     double distance = si_thickness*fabs(slope)/si_strip_pitch;
02779     distance = distance - 1.;
02780     int idistance = (int) distance;
02781     if (idistance < 1) idistance = 1;
02782 
02783     if (size < idistance) icluster = 0;
02784     else if (size == idistance) icluster = 1;
02785     else if (size > idistance) icluster = 2;
02786 
02787     if (icluster == 0 && size >=2 && idistance >=2) icluster = 1;
02788 
02789     return icluster;
02790 }
02791 
02792 //#########################################################################
02793 bool GFtutor::neighbourTowers(int itower, int jtower)
02794 //#########################################################################
02795 {
02796     bool ok = false;
02797 
02798     int kxtower = (int) itower/10;
02799     int kytower = itower - 10*itower;
02800 
02801     int kx0tower = (int) jtower/10;
02802     int ky0tower = jtower - 10*kx0tower;
02803 
02804     if (kxtower >= kx0tower-1 && kxtower <= kx0tower+1) ok = true;      
02805     if (kytower >= ky0tower-1 && kytower <= ky0tower+1) ok = true;
02806 
02807     // problem detected on the ID hit - not used for the momemt - 05/16/99 JAH
02808     if (itower < 10 || jtower < 10) ok =true;
02809 
02810     return ok;
02811 }