GFparticle.cpp Source File

00001 // $Id: GFparticle.cpp,v 1.4 2001/10/10 15:35:01 usher 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 
00016 
00017 #include "TkrRecon/GFparticle.h"
00018 //-----------------------------------------------------
00019 //
00020 //   GFtrack
00021 //
00022 //-----------------------------------------------------
00023 //######################################################
00024 GFtrack::GFtrack(SiCluster::view axis,
00025                  double sigmaCut,
00026                  double energy, 
00027                  int ist, 
00028                  const Ray& testRay,
00029                  bool run)
00030   : GFbase(sigmaCut,energy,ist,testRay),
00031     m_axis(axis), m_status(EMPTY),  
00032     m_qbest(-1e6), m_gaps(0), m_istGaps(0), m_lstLayer(0), m_noisyHits(0),
00033     m_istNoisyHits(0)
00034 //#######################################################
00035 {
00036     ini();
00037     if (run == true) {
00038         doit();
00039         fit();
00040         if (empty()) clear();
00041     }
00042 }
00043 
00044 //##########################################
00045 void GFtrack::flagAllHits(int iflag)
00046 //##########################################
00047 {
00048     for(unsigned int i=0; i<kplanelist.size(); i++) {
00049         GFtutor::_DATA->flagHit(m_axis, kplanelist[i].getIDHit(), iflag);
00050     }
00051 }
00052 
00053 //##########################################
00054 void GFtrack::unFlagAllHits()
00055 //##########################################
00056 {
00057     for (unsigned i=0; i<kplanelist.size(); i++) {
00058         GFtutor::_DATA->unflagHit(m_axis, kplanelist[i].getIDHit());
00059     }
00060 }  
00061 //########################################################
00062 bool GFtrack::empty() const
00063 //########################################################
00064 {
00065     bool empty = GFdata::empty();
00066     if (firstLayer() < 0) empty = true;
00067     if (numDataPoints() < GFcontrol::minSegmentHits) empty = true;
00068     if (chiSquare() < 0.) empty = true;
00069     return empty;
00070 }
00071 
00072 //########################################################
00073 bool GFtrack::accept() const
00074 //########################################################
00075 {
00076     bool valid = false;
00077     if (empty()) return valid;
00078     
00079     if (chiSquare() > GFcontrol::maxChiSq) return valid;
00080     if (Q() < GFcontrol::minQ) return valid;
00081     int idhit;
00082     double sigma;
00083     if (GFtutor::CUT_veto) {
00084         if (veto(idhit,sigma)) return valid;
00085     }
00086     
00087     valid = true;
00088     return valid;
00089 }
00090 //##########################################
00091 void GFtrack::clear()
00092 //##########################################
00093 {   
00094     
00095     KalTrack::clear();
00096     
00097     m_lstGaps      = 0;
00098     m_runChiSquare = 0.;
00099     m_gaps         = 0;
00100     m_istGaps      = 0;
00101     m_lstLayer     = 0;
00102     m_noisyHits    = 0;
00103     m_istNoisyHits = 0;
00104     
00105     m_status       = EMPTY;
00106 
00107     if (_mGFsegment) _mGFsegment->clear();
00108     
00109     m_qbest = -1e6;
00110     GFdata::ini();
00111     setAlive();
00112     
00113 }
00114 
00115 //########################################################
00116 void GFtrack::writeOut(MsgStream& log) const
00117 //########################################################
00118 {
00119     // kludge to avoid egcs warning messages
00120     int axis = getAxis();
00121     int status = m_status;
00122     log << MSG::DEBUG << " --- GFtrack::writeOut --- " << endreq;
00123     log << MSG::DEBUG << " axis           = " << axis << endreq;
00124     log << MSG::DEBUG << " Qbest          = " << Qbest() << endreq;
00125     log << MSG::DEBUG << " last  Layer    = " << lastLayer() << endreq;
00126     log << MSG::DEBUG << " num Hits       = " << numDataPoints() << endreq;
00127     log << MSG::DEBUG << " num Gaps       = " << numGaps() << endreq;
00128     log << MSG::DEBUG << " num First Gaps = " << numFirstGaps() << endreq;
00129     log << MSG::DEBUG << " num Noise      = " << numNoise() << endreq;
00130     log << MSG::DEBUG << " num First Noise= " << numFirstNoise() << endreq;
00131     log << MSG::DEBUG << " last Status    = " << status << endreq; 
00132     
00133 //    GFdata::writeOut(out);
00134     
00135 //    std::cout << " --> KalTrack : " << endreq;
00136 //    KalTrack::writeOut(out);
00137     
00138 }
00139 
00140 //########################################################
00141 void GFtrack::draw(gui::DisplayRep& v) 
00142 //########################################################
00143 {
00144         KalTrack::drawChiSq(v,m_axis,KalHit::SMOOTH);
00145         KalTrack::drawTrack(v,m_axis,KalHit::SMOOTH);
00146 }
00147 //#########################################################################
00148 bool GFtrack::veto(int& idhit, double& sigma) const
00149 //#########################################################################
00150 { 
00151     // WORK
00152     bool veto = false;
00153     
00154     int klayer = firstKPlane().getIDPlane() -1;
00155     if (klayer < 0) return veto;
00156     _mGFsegment->previous(klayer);
00157     
00158     if (_mGFsegment->status() == FOUND) {
00159         idhit = _mGFsegment->indexhit();
00160         sigma = _mGFsegment->getKPlane().getSigma(KalHit::PRED);
00161         if (sigma < GFcontrol::sigmaVeto) veto = true;
00162     }
00163     
00164     return veto;
00165 }
00166 
00167 //--------------------------------------------------------
00168 //  GFtrack - Private 
00169 //--------------------------------------------------------
00170 
00171 //########################################################
00172 void GFtrack::ini()
00173 //########################################################
00174 {
00175     m_status = EMPTY;
00176 
00177     m_gaps         = 0;
00178     m_istGaps      = 0;
00179     m_runChiSquare = 0.;
00180     m_lstGaps      = 0; 
00181     m_noisyHits    = 0;
00182     m_istNoisyHits = 0;
00183     m_lstLayer     = 0;
00184     
00185     KalTrack::clear();
00186     
00187     _mGFsegment = 0;
00188     _mGFsegment = new GFsegment(this);
00189     
00190     m_qbest = -1e6;
00191     GFdata::ini();
00192     setAlive();
00193     
00194 }
00195 
00196 //#######################################################
00197 void GFtrack::step(int kplane)
00198 //#######################################################
00199 {
00200     if (!alive()) return;
00201     
00202     m_status = EMPTY;
00203 
00204     if (kplane > GFtutor::numPlanes()-1) return;
00205     
00206     if (numDataPoints() == 0)
00207     {
00208         _mGFsegment->best(kplane);
00209     }
00210     else 
00211     {
00212         _mGFsegment->next(kplane);
00213     }
00214     
00215     m_status = _mGFsegment->status();
00216 
00217     if (m_status == FOUND) 
00218     {
00219         double segChiSquare = _mGFsegment->chiSquare();
00220         double maxChiSquare = numDataPoints() * 10. * m_runChiSquare;
00221 
00222         //If we have fewer than 4 hits on a track then we go ahead and add it
00223         //since, in theory, we can reject this segment later.
00224         //If the segment chi-square is small (< 1) then add the hit as well (the
00225         //chi-square does not appear to be well determined because the energy is not
00226         //well determined so for stiff tracks the chi-square can be effectively zero)
00227         //Finally, add the hit if the chi-square doesn't blow up too badly according to 
00228         //a sliding scale which depends on track length.
00229         //I dearly hope the "new" tracking code does a better job at this!
00230        if (numDataPoints() < 4 || segChiSquare < 1. || segChiSquare < maxChiSquare)
00231         {
00232             kplanelist.push_back(_mGFsegment->getKPlane());
00233             if (numDataPoints() > 2) m_runChiSquare += segChiSquare;
00234         }
00235         else 
00236         {
00237             setStatus(EMPTY);
00238         }
00239     }
00240     
00241 }
00242 
00243 //#########################################################
00244 void GFtrack::anastep(int kplane) 
00245 //#########################################################
00246 {
00247     if (!alive()) return;
00248 
00249     contability(kplane);
00250     if (end()) {
00251         kill();
00252     }
00253 }
00254 
00255 //##########################################################
00256 void GFtrack::contability(int kplane) 
00257 //##########################################################
00258 {
00259     if (!alive()) return;
00260     
00261     if (m_status != FOUND) {
00262         m_lstGaps++;
00263         if (numDataPoints()>0 && numDataPoints()<3) m_istGaps++;
00264     }
00265     
00266     // if (m_statushit == CRACK) m_cracks++;
00267     if (m_status == EMPTY) m_gaps++;
00268     
00269     if (m_status == FOUND && numDataPoints() > 0) {
00270         m_lstGaps =0;
00271 
00272         m_lstLayer = lastKPlane().getIDPlane();
00273         int indexHit = lastKPlane().getIDHit();
00274         /*int type = GFtutor::_DATA->clusterNoise(m_axis, indexHit);
00275         if(type > 0) {
00276             m_noisyHits++;
00277             if (kplanelist.size() < 3 ) m_noisyHits++;
00278         }*/
00279     } 
00280     
00281    // if (m_status == FOUND && numDataPoints() == 0) CONTROL_error++;
00282     
00283 }
00284 
00285 //#########################################################################
00286 void GFtrack::fit()
00287 //#########################################################################
00288 {
00289     
00290     GFdata::ini();
00291     
00292     if (kplanelist.size()< 3) {
00293         KalTrack::clear();
00294         return;
00295     }
00296     
00297     //----------------------------
00298     // Voila monsieur Kalman
00299     //----------------------------
00300     KalTrack::doFit();
00301     //----------------------------
00302     
00303     loadGFdata();
00304 }
00305 //#########################################################################
00306 void GFtrack::loadGFdata()
00307 //#########################################################################
00308 {
00309     
00310     // m_qbest = doQbest();
00311     
00312     // Output Data
00313     double x0,y0,z0;
00314     double slopex,slopey;
00315     x0=y0=z0=0.;
00316     slopex=slopey=0.;
00317     z0=kplanelist[0].getZPlane();
00318     if (m_axis == SiCluster::X) {
00319         x0=position(0.);
00320         slopex=slope();
00321     } else {
00322         y0=position(0.);
00323         slopey=slope();
00324     }
00325     m_vertex = Point(x0,y0,z0);
00326     double factor = -1.;
00327     m_direction = Vector(factor*slopex,factor*slopey,factor).unit();
00328     m_RCenergy =  KalEnergy();
00329     m_quality = computeQuality();
00330     m_firstLayer = kplanelist[0].getIDPlane();
00331     m_nhits = numDataPoints();
00332     m_itower = kplanelist[0].getIDTower();
00333     
00334 } 
00335 
00336 //#########################################################################
00337 double GFtrack::computeQuality() const
00338 //#########################################################################
00339 {
00340     double quality = 0.;
00341     quality = 27./(9.+chiSquare()) + 15./(5.+chiSquareSegment()) + 
00342         2.*(numDataPoints()-3.) - numGaps() - 5.*numFirstGaps();
00343     //          + (1./(0.1+abs(kink(0))));
00344     return quality;
00345 }
00346 //#########################################################################
00347 bool GFtrack::end() const
00348 //#########################################################################
00349 {
00350     bool end = !alive();
00351     if (m_lstGaps >= GFcontrol::maxConsecutiveGaps) end = true;
00352     return end;
00353 }                 
00354 
00355 //#########################################################################
00356 void GFtrack::kill()
00357 //#########################################################################
00358 {
00359     if (m_status == EMPTY) {
00360         m_gaps+=(-m_lstGaps);
00361         m_lstGaps = 0;
00362     }
00363     m_status = EMPTY;
00364 
00365     m_alive = false;
00366     
00367 }
00368 //########################################################
00369 void GFtrack::setAlive()
00370 //########################################################
00371 {
00372     m_alive = true;
00373 }
00374 //########################################################
00375 void GFtrack::setIniEnergy(double ene)
00376 //########################################################
00377 {
00378     m_iniEnergy = ene;
00379     KalTrack::setIniEnergy(ene);
00380 }
00381 
00382 //################################################
00383 KalPlane GFtrack::firstKPlane() const
00384 //################################################
00385 {
00386     if (kplanelist.size() == 0) {
00387         std::cout << "ERROR GFtrack::thisKPlane " << endreq;
00388         return originalKPlane();
00389     }
00390     return kplanelist.front();
00391 }
00392 
00393 //################################################
00394 KalPlane GFtrack::lastKPlane() const
00395 //################################################
00396 {
00397     if (kplanelist.size() == 0) {
00398         return originalKPlane();
00399     }
00400     return kplanelist.back();
00401 }
00402 //################################################
00403 KalPlane GFtrack::previousKPlane() const 
00404 //################################################
00405 {
00406     if (kplanelist.size() <= 1) {
00407         //              std::cout << "ERROR GFtrack::previousKPlane " << endreq;
00408         return originalKPlane();
00409     }
00410     int iprevious = kplanelist.size()-2;
00411     if (iprevious == -1) return originalKPlane();
00412     return kplanelist[iprevious];
00413 }
00414 //################################################
00415 KalPlane GFtrack::originalKPlane() const
00416 //################################################
00417 {
00418     
00419     Ray testRay(m_inVertex,m_inDirection);
00420     
00421     double x0=testRay.position(0.).x();
00422     double y0=testRay.position(0.).y();
00423     double z0=testRay.position(0.).z();
00424     
00425     double dirX=testRay.direction().x();
00426     double dirY=testRay.direction().y();
00427     double dirZ=testRay.direction().z();
00428     
00429     double x = 0.;
00430     double x_orth = 0.;
00431     m_axis == SiCluster::X? x = x0: x=y0;
00432     m_axis == SiCluster::X? x_orth = y0: x_orth = x0;
00433     
00434     double slope=0.;
00435     double slope_orth=0.;
00436     if (dirZ != 0.) {
00437         slope = (m_axis == SiCluster::X? dirX/dirZ: dirY/dirZ);
00438         slope_orth = (m_axis == SiCluster::X? dirY/dirZ: dirX/dirZ);
00439     }
00440     
00441     KalPar porth(x_orth, slope_orth);
00442     
00443     KalPar pfit(x,slope);
00444     
00445     //unused:    double sigma = SiTracker::siStripPitch()/sqrt(12.);
00446     /*KalMatrix Q=KalPlane::Q(m_iniEnergy,slope,slope_orth, 
00447     KalPlane::radLen(m_iniLayer)); */
00448     double sigma2Slope = GFcontrol::iniErrorSlope * GFcontrol::iniErrorSlope;
00449     double sigma2Position = GFcontrol::iniErrorPosition * GFcontrol::iniErrorPosition;
00450     KalMatrix covfit(sigma2Position,sigma2Slope,0.);
00451     
00452     KalHit hitfit(KalHit::FIT, pfit, covfit);
00453     
00454     KalPlane kp(-1,m_iniLayer-1,m_iniEnergy, z0, porth, hitfit);
00455     
00456     return kp;
00457 }
00458 
00459 //#######################################################
00460 void GFtrack::removeStep(int iplane)
00461 //#######################################################
00462 {
00463     if (iplane == -1) iplane = numDataPoints()-1;
00464  //   if (iplane == -1) CONTROL_error++;
00465     
00466     if (iplane == numDataPoints()-1) {
00467         kplanelist.pop_back();
00468     } else {
00469         // WORK : remove the plane k 
00470     }
00471     
00472     setStatus(EMPTY);
00473 }
00474 
00475 //#########################################################
00476 double GFtrack::doQbest()
00477 //#########################################################
00478 {
00479     double qbest = -1e6;
00480     if (numDataPoints()<3) return qbest;
00481     KalTrack::doFit();
00482     loadGFdata();
00483     // qbest=-1.*chiSquareSegment(m_sigmaCut*m_sigmaCut);
00484     m_qbest = Q();
00485     GFdata::ini();
00486     return m_qbest;
00487 }
00488 
00489 //#########################################################################
00490 void GFtrack::associateOrthStep(const GFtrack* _GFtrk, KalHit::TYPE typ)
00491 //#########################################################################
00492 {
00493     // locate the FIT hit of the previous plane in this step!
00494     if (numDataPoints() >=1 && _GFtrk->numDataPoints() >=1 && m_status == FOUND) {
00495 
00496         KalPar XPar = _GFtrk->lastKPlane().getHit(typ).getPar();
00497 
00498         double dz = kplanelist.back().getZPlane()-
00499             _GFtrk->kplanelist.back().getZPlane();
00500         double x0 = XPar.getPosition()+dz*XPar.getSlope();
00501         KalPar OrthPar(x0,XPar.getSlope());
00502 
00503         kplanelist.back().setOrthPar(OrthPar); // it changes the plane contents
00504     }
00505 }
00506 
00507 //#########################################################################
00508 void GFtrack::associateOrthGFtrack(const GFtrack* _GFtrk, bool fix, KalHit::TYPE typ)
00509 //#########################################################################
00510 {
00511     if (numDataPoints() == 0 || _GFtrk->numDataPoints() == 0) return;
00512     // Associate another track
00513     double XLIMIT = -91.;
00514     double SLOPENULL = 0.;
00515     int mplanes = _GFtrk->numDataPoints();
00516     for ( int iplane = 0; iplane < numDataPoints(); iplane++) {
00517         int idplane = kplanelist[iplane].getIDPlane();
00518         int jplane = -1;
00519         int jdplane = 0;
00520         do {
00521             jplane++;
00522             jdplane = _GFtrk->kplanelist[jplane].getIDPlane();
00523         } while (jdplane <= idplane && jplane+1 < mplanes-1);
00524         if (jdplane > idplane && jplane > 0) jplane--;
00525         
00526         
00527         KalPar Ypar = _GFtrk->kplanelist[jplane].getHit(typ).getPar();
00528         double dz = kplanelist[iplane].getZPlane()
00529             - _GFtrk->kplanelist[jplane].getZPlane();
00530         double x0 = Ypar.getPosition()+dz*Ypar.getSlope();
00531         double slope0 = Ypar.getSlope();
00532         if (!fix) {
00533             x0 = XLIMIT;
00534             slope0 = SLOPENULL;
00535         }
00536         KalPar OrthPar(x0,slope0);
00537 
00538         kplanelist[iplane].setOrthPar(OrthPar);
00539     }
00540 }
00541 //-----------------------------------------------------
00542 //
00543 //   GFparticle
00544 //
00545 //-----------------------------------------------------
00546 //######################################################
00547 GFparticle::GFparticle(double sigmaCut,
00548                        double energy, 
00549                        int ist, 
00550                        const Ray& testRay,
00551                        bool run) : GFbase(sigmaCut,energy,ist,testRay)
00552                        //#######################################################
00553 {
00554     ini();
00555     if (run == true) {
00556         doit();
00557         fit();
00558         if (empty()) clear();
00559     }
00560 }
00561 
00562 //##########################################
00563 void GFparticle::flagAllHits(int iflag)
00564 //##########################################
00565 {
00566     _mXGFtrack->flagAllHits(iflag);
00567     _mYGFtrack->flagAllHits(iflag);
00568 }
00569 
00570 //##########################################
00571 void GFparticle::unFlagAllHits()
00572 //##########################################
00573 {
00574     _mXGFtrack->unFlagAllHits();
00575     _mYGFtrack->unFlagAllHits();
00576 }  
00577 //########################################################
00578 bool GFparticle::empty() const
00579 //########################################################
00580 {
00581     bool empty = GFdata::empty();
00582     if (empty) return empty;
00583     empty = _mXGFtrack->empty() || _mYGFtrack->empty();
00584     return empty;
00585 }
00586 
00587 //########################################################
00588 bool GFparticle::accept() const
00589 //########################################################
00590 {
00591     bool valid = false;
00592     if (empty()) return valid;
00593     
00594     if (Q() < GFcontrol::minQ) return valid;
00595     
00596     if (GFtutor::CUT_veto) {
00597         int idhitX = -1;
00598         double sigmaX = -1.;
00599         int idhitY = -1;
00600         double sigmaY = -1.;
00601         if (_mXGFtrack->veto(idhitX,sigmaX) && 
00602             _mYGFtrack->veto(idhitY,sigmaY)) return valid;
00603     }
00604     
00605     valid = true;
00606     return valid;
00607 }
00608 //##########################################
00609 void GFparticle::clear()
00610 //##########################################
00611 {   
00612     _mXGFtrack->clear();
00613     _mYGFtrack->clear();
00614     
00615     m_gaps   = 0;
00616     m_istGaps= 0;
00617     m_noisyHits = 0;
00618     m_istNoisyHits   = 0;
00619     m_lstLayer = 0;
00620     
00621     m_status = EMPTY;
00622     m_associate = true;
00623     m_conflictPattern = false;
00624     
00625     m_qbest = -1e6;
00626     GFdata::ini();
00627     setAlive();
00628     
00629 }
00630 
00631 //########################################################
00632 void GFparticle::writeOut(MsgStream& log) const
00633 //########################################################
00634 {
00635     // kludge to avoid egcs warning messages
00636     int status = m_status;
00637     log << MSG::DEBUG << " --- GFparticle::writeOut --- " << endreq;
00638     log << MSG::DEBUG << " Qbest          = " << Qbest() << endreq;
00639     log << MSG::DEBUG << " last  Layer    = " << lastLayer() << endreq;
00640     log << MSG::DEBUG << " num Gaps       = " << numGaps() << endreq;
00641     log << MSG::DEBUG << " num First Gaps = " << numFirstGaps() << endreq;
00642     log << MSG::DEBUG << " num Noise      = " << numNoise() << endreq;
00643     log << MSG::DEBUG << " num First Noise= " << numFirstNoise() << endreq;
00644     log << MSG::DEBUG << " last Status    = " << status << endreq; 
00645     
00646     GFdata::writeOut(log);
00647     
00648     log << MSG::DEBUG << " -->  X - GFtrack : " << endreq;
00649     _mXGFtrack->writeOut(log);
00650     log << MSG::DEBUG << " -->  Y - GFtrack : " << endreq;
00651     _mYGFtrack->writeOut(log);
00652 }
00653 
00654 //########################################################
00655 void GFparticle::draw(gui::DisplayRep& v) 
00656 //########################################################
00657 {       
00658         v.setColor("black");
00659         _mXGFtrack->draw(v);
00660         _mYGFtrack->draw(v);
00661 }
00662 //--------------------------------------------------------
00663 //  GFparticle - Private 
00664 //--------------------------------------------------------
00665 
00666 //########################################################
00667 void GFparticle::ini()
00668 //########################################################
00669 {
00670     m_gaps   = 0;
00671     m_istGaps= 0;
00672     m_noisyHits = 0;
00673     m_istNoisyHits   = 0;
00674     m_lstLayer = 0;
00675     
00676     m_status = EMPTY;
00677     m_associate = true;
00678     m_conflictPattern = false;
00679     
00680     Ray testRay(m_inVertex,m_inDirection);      
00681     _mXGFtrack = new GFtrack(SiCluster::X, m_sigmaCut, 
00682         m_iniEnergy, m_iniLayer, testRay, false);
00683     
00684     _mYGFtrack = new GFtrack(SiCluster::Y, m_sigmaCut,
00685         m_iniEnergy, m_iniLayer, testRay, false);
00686     
00687     m_qbest = -1e6;
00688     GFdata::ini();
00689     setAlive();
00690     
00691 }
00692 
00693 //#######################################################
00694 void GFparticle::step(int kplane)
00695 //#######################################################
00696 {
00697     
00698     if (!alive()) return;
00699     
00700     _mXGFtrack->step(kplane);
00701     _mYGFtrack->step(kplane);
00702     
00703     if (m_associate) {          
00704         associateStatus();
00705         associateStep();
00706     }
00707 }
00708 //#########################################################################
00709 void GFparticle::anastep(int kplane) 
00710 //#########################################################################
00711 {
00712     if (!alive()) return;
00713     
00714     _mXGFtrack->anastep(kplane);
00715     _mYGFtrack->anastep(kplane);
00716     
00717     contability(kplane);
00718     
00719     if (m_associate) {
00720         associateAnaStep();
00721     }
00722     
00723     if (end()) {
00724         kill();
00725     }
00726 }
00727 
00728 //#########################################################################
00729 void GFparticle::fit()
00730 //#########################################################################
00731 {
00732     GFdata::ini();
00733     
00734     _mXGFtrack->fit();
00735     _mYGFtrack->fit();
00736     
00737     if (!_mXGFtrack->empty() && !_mYGFtrack->empty()) 
00738         loadGFdata();
00739     
00740     if (m_associate && 
00741         (!_mXGFtrack->empty() && !_mYGFtrack->empty())) 
00742         associateFit();
00743 }
00744 
00745 //#########################################################################
00746 bool GFparticle::end() const
00747 //#########################################################################
00748 {
00749     bool end = !alive();
00750     // This comparatios is between different status and makes no sense
00751     // if (m_status == DONE) return end = true;
00752     if (!_mXGFtrack->alive() || !_mYGFtrack->alive()) end = true;
00753     return end;
00754 }
00755 
00756 //#########################################################################
00757 void GFparticle::kill()
00758 //#########################################################################
00759 {
00760     m_alive = false;
00761     _mXGFtrack->kill();
00762     _mYGFtrack->kill();
00763     
00764 }
00765 //#########################################################################
00766 void GFparticle::setAlive()
00767 //#########################################################################
00768 {
00769     m_alive = true;
00770     _mXGFtrack->setAlive();
00771     _mYGFtrack->setAlive();
00772     
00773 }
00774 
00775 
00776 //#########################################################################
00777 void GFparticle::loadGFdata()
00778 //#########################################################################
00779 {
00780     int ixlayer = _mXGFtrack->firstLayer();
00781     int iylayer = _mYGFtrack->firstLayer();
00782     
00783     m_firstLayer = (ixlayer <= iylayer? ixlayer : iylayer);
00784     
00785     int ixtower = _mXGFtrack->tower();
00786     //    int iytower = _mYGFtrack->tower();
00787     m_itower = ixtower;
00788     
00789     m_nhits = _mXGFtrack->nhits() + _mYGFtrack->nhits();
00790     
00791     m_quality = _mXGFtrack->Q() + _mYGFtrack->Q();
00792     
00793     m_RCenergy = 0.5*(_mXGFtrack->RCenergy() + _mYGFtrack->RCenergy());
00794     
00795     Ray XRay = _mXGFtrack->ray();
00796     Ray YRay = _mYGFtrack->ray();
00797     m_vertex=GFbase::doVertex(XRay, YRay);
00798     
00799     m_direction = GFbase::doDirection(_mXGFtrack->direction(),_mYGFtrack->direction());  
00800 }
00801 //#########################################################################
00802 void GFparticle::contability(int kplane) 
00803 //#########################################################################
00804 {
00805     if (!alive()) return;
00806 }
00807 
00808 //#########################################################################
00809 void GFparticle::associateStep() 
00810 //#########################################################################
00811 {
00812     bool ok = false;
00813     bool done = false;
00814 
00815     ok = GFparticle::sameTower(_mXGFtrack,_mYGFtrack);
00816     if (!ok) done = GFparticle::removeWorseStep(_mXGFtrack,_mYGFtrack);
00817     
00818     if (!ok && !done) m_conflictPattern = true;
00819     
00820 }
00821 
00822 //#########################################################################
00823 void GFparticle::associateStatus() 
00824 //#########################################################################
00825 {
00826     if (_mXGFtrack->status() == FOUND || _mYGFtrack->status() == FOUND)
00827         m_status = FOUND;
00828     else m_status = EMPTY;
00829 }
00830 
00831 //#########################################################################
00832 void GFparticle::associateAnaStep() 
00833 //#########################################################################
00834 {
00835     _mXGFtrack->associateOrthStep(_mYGFtrack);  
00836     _mYGFtrack->associateOrthStep(_mXGFtrack);
00837 }
00838 
00839 //#########################################################################
00840 void GFparticle::associateFit() 
00841 //#########################################################################
00842 {
00843     _mXGFtrack->associateOrthGFtrack(_mYGFtrack, m_associate, KalHit::SMOOTH);  
00844     _mYGFtrack->associateOrthGFtrack(_mXGFtrack, m_associate, KalHit::SMOOTH);
00845 }
00846 
00847 //#########################################################################
00848 bool GFparticle::sameTower(const GFtrack* _GFtrack1, const GFtrack* _GFtrack2) 
00849 //#########################################################################
00850 {
00851     bool sametower = true;
00852     int itower1 = 0;
00853 
00854     if (_GFtrack1->status() == FOUND && _GFtrack1->numDataPoints() > 0) {
00855         itower1 = _GFtrack1->lastKPlane().getIDTower();
00856     }
00857     int itower2 = 0;
00858 
00859     if (_GFtrack2->status() == FOUND && _GFtrack2->numDataPoints() > 0) {
00860         itower2 = _GFtrack2->lastKPlane().getIDTower();
00861     }
00862     
00863     if (itower1 >= 10 && itower2 >= 10 && itower1 != itower2) sametower = false;
00864     return sametower;
00865 }
00866 //#########################################################################
00867 bool GFparticle::removeWorseStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2) 
00868 //#########################################################################
00869 {
00870     bool done = false;
00871     
00872 //    if (_GFtrack1->status() != FOUND || _GFtrack2->status() != FOUND) CONTROL_error++;
00873 //    if (_GFtrack1->numDataPoints() == 0 || _GFtrack2->numDataPoints() == 0) CONTROL_error++;
00874     
00875     done = true;
00876     
00877     int kplane1 = _GFtrack1->previousKPlane().getIDPlane();
00878     int kplane2 = _GFtrack2->previousKPlane().getIDPlane();
00879     
00880     // the closest track first
00881     if (kplane1 != kplane2) {
00882         if (kplane1 > kplane2) _GFtrack2->removeStep();
00883         else  _GFtrack1->removeStep();
00884         return done;
00885     }
00886     
00887     int nhits1 = _GFtrack1->numDataPoints();
00888     int nhits2 = _GFtrack2->numDataPoints();
00889     
00890     double chi1 = _GFtrack1->lastKPlane().getDeltaChiSq(KalHit::FIT);
00891     double chi2 = _GFtrack2->lastKPlane().getDeltaChiSq(KalHit::FIT);
00892     if (nhits1 < 2 || nhits2 < 2) {
00893         if (_GFtrack1->_mGFsegment->accept() || 
00894             _GFtrack2->_mGFsegment->accept()  ) {
00895             if (_GFtrack1->_mGFsegment->accept()) chi1 = _GFtrack1->_mGFsegment->chiGFSq();
00896             else chi1 = 1e6;
00897             if (_GFtrack2->_mGFsegment->accept()) chi2 = _GFtrack2->_mGFsegment->chiGFSq();
00898             else chi2 = 1e6;
00899         }
00900     } 
00901     
00902     if (chi2 > chi1) _GFtrack2->removeStep();
00903     else if (chi2 < chi1 ) _GFtrack1->removeStep();
00904     else done = false;
00905     
00906     return done;
00907 }