GFgamma.cpp Source File

00001 //-------------------------------------------------------------------------
00002 //
00003 //     GFpair 
00004 //
00005 //-------------------------------------------------------------------------
00006 
00007 #include "TkrRecon/GFgamma.h"
00008 
00009 //#########################################################################
00010 GFpair::GFpair(double xene, enum SiCluster::view axis,
00011                double sigmaCut,
00012                double energy, 
00013                int ist, 
00014                const Ray& testRay,
00015                bool run) : GFbase(sigmaCut,energy,ist,testRay),
00016                m_xEne(xene),
00017                m_axis(axis)
00018 //##########################################################################
00019 {
00020     ini();
00021     if (run == true) {
00022                 doit();
00023                 fit();
00024                 if (empty()) clear();
00025     }
00026 }
00027 //#########################################################################
00028 void GFpair::flagAllHits(int iflag)
00029 //#########################################################################
00030 {
00031     _mGFbest->flagAllHits(iflag);
00032     _mGFpair->flagAllHits(iflag);
00033 }
00034 //#########################################################################
00035 void GFpair::unFlagAllHits()
00036 //#########################################################################
00037 {
00038     _mGFbest->unFlagAllHits();
00039     _mGFpair->unFlagAllHits();
00040 }
00041 //#########################################################################
00042 bool GFpair::empty() const
00043 //#########################################################################
00044 {
00045     bool empty = GFdata::empty();
00046     if (empty) return empty;
00047     
00048     empty = empty || _mGFbest->empty();
00049     return empty;
00050 }
00051 
00052 //#########################################################################
00053 bool GFpair::accept() const
00054 //#########################################################################
00055 {
00056     bool ok = false;
00057     if (_mGFbest->empty()) return ok;
00058     
00059     // if there is a pair please do not impose the veto!
00060     int indexhit;
00061     double sigma;
00062     ok = true;
00063     if (GFtutor::CUT_veto) ok = !_mGFbest->veto(indexhit,sigma);
00064     
00065     return ok;
00066 }
00067 
00068 //#########################################################################
00069 void GFpair::clear()
00070 //#########################################################################
00071 {
00072     m_status = TOGETHER;
00073     m_decideBest = false;
00074     
00075     m_weightBest = 0.;
00076     m_errorSlope = 0.;
00077     
00078     m_together = 0;
00079     m_split = 0;
00080     m_one = 0;
00081     m_shared = 0;
00082     m_empty = 0;
00083     
00084     _mGFbest->clear();
00085     _mGFpair->clear();  
00086     if (_mGFalive) _mGFalive->clear();
00087     
00088     GFdata::ini();
00089     setAlive();
00090 }
00091 
00092 //########################################################
00093 void GFpair::writeOut(MsgStream& log) const
00094 //########################################################
00095 {
00096     // kludge to avoid warnings from egcs
00097     int axis   = m_axis;
00098     int status = m_status;
00099     log << MSG::DEBUG << " --- GFpair::writeOut --- " << "\n";
00100     log << MSG::DEBUG << " axis:           " << axis << "\n";
00101     log << MSG::DEBUG << " Energy Split    " << m_xEne << "\n";
00102     log << MSG::DEBUG << " Weight          " << m_weightBest << "\n";
00103     log << MSG::DEBUG << " planes together " << numTogether() << "\n";
00104     log << MSG::DEBUG << " planes split    " << numSplit() << "\n";
00105     log << MSG::DEBUG << " planes one      " << numOne() << "\n";
00106     log << MSG::DEBUG << " shared hits     " << numSharedHits() << "\n";
00107     log << MSG::DEBUG << " empty planes    " << numEmpty() << "\n";
00108     log << MSG::DEBUG << " last Status     " << status << "\n"; 
00109     
00110     GFdata::writeOut(log);
00111 //    log << MSG::DEBUG << " --> best Track : " << "\n";
00112 //    _mGFbest->writeOut(out);
00113 //    log << MSG::DEBUG << " --> pair Track : " << "\n";
00114 //    _mGFpair->writeOut(out);
00115 }
00116 
00117 //########################################################
00118 void GFpair::draw(gui::DisplayRep& v) 
00119 //########################################################
00120 {
00121         v.setColor("blue");
00122         _mGFbest->draw(v);
00123         v.setColor("green");
00124         _mGFpair->draw(v);
00125 }
00126 //-------------------------------------------------------------------------
00127 //    GFpair - Private
00128 //-------------------------------------------------------------------------
00129 //#########################################################################
00130 void GFpair::ini()
00131 //#########################################################################
00132 {
00133     // status
00134     m_status = TOGETHER;
00135     m_decideBest = false;
00136     
00137     // contability
00138     m_together = 0;
00139     m_split = 0;
00140     m_one = 0;
00141     m_shared = 0;
00142     m_empty = 0;
00143     
00144     // pair addition results
00145     m_weightBest = 0.;
00146     m_errorSlope = 0.;
00147     
00148     // internal variables
00149     _mGFbest = 0;
00150     _mGFpair = 0;
00151     _mGFalive = 0;
00152     
00153     Ray testRay(m_inVertex,m_inDirection);      
00154 
00155     _mGFbest = new GFtrack(m_axis, sigmaCut(), 
00156                 GFcontrol::FEne*m_iniEnergy, m_iniLayer, testRay, false);
00157     
00158     _mGFpair = new GFtrack(m_axis,  sigmaCut(),
00159                 (1.-GFcontrol::FEne)*m_iniEnergy, m_iniLayer, testRay, false);
00160     
00161     GFdata::ini();
00162     setAlive();
00163 }
00164 
00165 //#########################################################################
00166 void GFpair::step(int kplane)
00167 //#########################################################################
00168 {    
00169     if (!alive()) return;
00170     
00171     newStatus(kplane);
00172     
00173     switch (status()) {
00174     case TOGETHER:
00175         stepTogether(kplane);
00176         break;
00177     case SPLIT:
00178         stepSplit(kplane);
00179         break;
00180     case ONE:
00181         _mGFalive->step(kplane);
00182         break;
00183     case DONE:
00184         break;
00185     }
00186 }
00187 
00188 //#########################################################################
00189 void GFpair::anastep(int kplane)
00190 //#########################################################################
00191 {       
00192     if (!m_alive) return;
00193 
00194     _mGFbest->anastep(kplane);
00195     _mGFpair->anastep(kplane);
00196     
00197     contability(kplane);
00198     
00199     if (end()) {
00200         kill();
00201     }
00202 }
00203 //#########################################################################
00204 void GFpair::fit()
00205 //#########################################################################
00206 {
00207     GFdata::ini();
00208     // Decide wich is the best track
00209     if (!m_decideBest) decideBest();
00210     
00211     if (m_decideBest) setIniEnergy();
00212     _mGFbest->fit();
00213     _mGFpair->fit();
00214     
00215     if (!_mGFbest->empty()) loadGFdata();
00216 }
00217 
00218 //#########################################################################
00219 bool GFpair::end() const
00220 //#########################################################################
00221 {
00222     bool end = !alive();
00223     if (m_status == DONE) return end = true;
00224     if (!_mGFbest->alive() && !_mGFpair->alive()) end = true;
00225     return end;
00226 }
00227 //#########################################################################
00228 void GFpair::kill() 
00229 //#########################################################################
00230 {
00231     m_alive = false;
00232     _mGFbest->kill();
00233     _mGFpair->kill();
00234 }
00235 //#########################################################################
00236 void GFpair::setAlive() 
00237 //#########################################################################
00238 {
00239     m_alive = true;
00240     _mGFbest->setAlive();
00241     _mGFpair->setAlive();
00242 }
00243 
00244 //#########################################################################
00245 void GFpair::contability(int kplane)
00246 //#########################################################################
00247 {
00248     if (m_status == TOGETHER) m_together++;
00249     else if (m_status == SPLIT) m_split++;
00250     else if (m_status == ONE) m_one++; 
00251     
00252     if (m_status == SPLIT) {
00253         if (_mGFbest->status() == FOUND && _mGFpair->status() == FOUND) {
00254 //          if (_mGFbest->numDataPoints() == 0 || _mGFpair->numDataPoints() == 0) GFcontrol::error++;
00255 //          else {
00256                 int idhit1 = _mGFbest->lastKPlane().getIDHit();
00257                 int idhit2 = _mGFpair->lastKPlane().getIDHit();
00258                 if (idhit1 == idhit2) m_shared++;
00259 //          }
00260         }
00261     } 
00262     
00263     if (m_status == TOGETHER) {
00264         if (_mGFbest->status() != FOUND) m_empty++;
00265     } else if (m_status == ONE) {
00266         if (_mGFalive->status() != FOUND) m_empty++;
00267     } else if (m_status == SPLIT) {
00268         if (_mGFbest->status() != FOUND && _mGFpair->status() != FOUND) m_empty++;
00269     } 
00270     
00271     if (m_status == ONE && m_split == 0 && m_together ==0) {
00272         m_one += m_together;
00273         m_together = 0;
00274     } 
00275     
00276 }
00277 //#########################################################################
00278 void GFpair::loadGFdata()
00279 //#########################################################################
00280 {
00281     m_firstLayer = _mGFbest->firstLayer();
00282     m_itower = _mGFbest->tower();
00283     m_nhits = _mGFbest->nhits();
00284     if (!_mGFpair->empty()) m_nhits += _mGFpair->nhits()-m_shared;
00285     
00286     if (!_mGFpair->empty()) {   
00287         m_quality = _mGFbest->Q() + _mGFpair->Q() 
00288             - (m_together/_mGFbest->numDataPoints())*(_mGFbest->Q());
00289     } else m_quality = _mGFbest->Q();
00290     
00291     m_RCenergy = doEnergy(_mGFbest, _mGFpair);
00292     
00293     if (_mGFpair->empty()) {
00294         m_direction = _mGFbest->direction();
00295         m_weightBest = 1.;
00296         m_errorSlope     = _mGFbest->errorSlopeAtVertex();
00297     } else {
00298                 if (GFcontrol::addTracksChi2) m_direction = doDirection(m_weightBest);
00299         else  m_direction = doDirection(_mGFbest, _mGFpair, m_weightBest, m_errorSlope);
00300         // m_direction = doDirectionXene(m_RCenergy, m_weightBest);
00301     }
00302     
00303     
00304     m_vertex = _mGFbest->vertex();
00305     if (!_mGFpair->empty()) {    
00306         Ray bestRay = Ray(_mGFbest->vertex(),_mGFbest->direction());
00307         Ray pairRay = Ray(_mGFpair->vertex(),_mGFpair->direction());
00308         m_vertex=GFbase::doVertex(bestRay, pairRay);
00309     }
00310 }
00311 //########################################################################
00312 void GFpair::newStatus(int kplane)
00313 //#########################################################################
00314 {
00315     StatusPair newStatus = status();
00316     
00317     switch(status()){
00318     case TOGETHER:
00319         if (!_mGFbest->alive()) newStatus = DONE;
00320         else if (forceSplit(kplane)) newStatus = ONE;
00321         break;
00322     case SPLIT:
00323         // You can be split but with only one hit, the together step should be used.
00324         if (_mGFbest->numDataPoints() <= 0) newStatus = TOGETHER; 
00325         if (!_mGFbest->alive() && !_mGFpair->alive()) newStatus = DONE;
00326         else if (!_mGFbest->alive() || !_mGFpair->alive()) newStatus = ONE;
00327         break;
00328     case ONE:
00329         if (!_mGFalive->alive()) newStatus = DONE;
00330         break;
00331     case DONE:
00332         break;
00333     }
00334     
00335     setStatus(newStatus);
00336 }
00337 //########################################################################
00338 bool GFpair::forceSplit(int kplane) const
00339 //########################################################################
00340 {
00341     // WORK : define the best spliting point according with the energy
00342     bool split = false;
00343     if ( kplane > GFcontrol::minSegmentHits + m_iniLayer -1) split = true;
00344     return split;
00345 }
00346 //########################################################################
00347 void GFpair::setStatus(StatusPair newStatus)
00348 //########################################################################
00349 {
00350     if (newStatus == status()) return;
00351     
00352     if (newStatus == ONE) {
00353         if (m_status == TOGETHER) {
00354             _mGFpair->clear();
00355             _mGFpair->kill();
00356         }
00357         _mGFalive = (_mGFbest->alive()? _mGFbest : _mGFpair);
00358     } else if (newStatus == DONE) {
00359         kill();
00360     }
00361     
00362     if (m_status == DONE || newStatus  == DONE) m_status = DONE;
00363     else if (m_status == ONE || newStatus == ONE) m_status = ONE;
00364     else if (m_status == SPLIT || newStatus == SPLIT) {
00365         m_status = SPLIT;
00366         if (_mGFbest->numDataPoints() == 0) m_status = TOGETHER;
00367     }
00368     else m_status = TOGETHER;
00369     
00370     // error control
00371 /*    if (m_status == ONE) {
00372         if (_mGFbest->alive() && _mGFpair->alive()) GFcontrol::error++;
00373     } else if (m_status != DONE) {
00374         if (!_mGFbest->alive() || !_mGFpair->alive()) GFcontrol::error++;
00375     } else {
00376         if (_mGFbest->alive() || _mGFpair->alive()) GFcontrol::error++;
00377     } */
00378 }
00379 //#########################################################################
00380 void GFpair::stepTogether(int klayer)
00381 //#########################################################################
00382 {
00383     _mGFbest->step(klayer);
00384     if (!_mGFbest->_mGFsegment->accept()) return;
00385     
00386     _mGFbest->_mGFsegment->flagUsedHits(klayer);
00387     _mGFpair->_mGFsegment->best(klayer);
00388     _mGFbest->_mGFsegment->unFlagUsedHits(klayer);
00389     
00390     bool split = _mGFpair->_mGFsegment->accept();
00391     if (split && _mGFpair->numDataPoints() == 0 ) {
00392         if (_mGFpair->_mGFsegment->getKPlane().getSigma(KalHit::SMOOTH) > sigmaCut()) {
00393             split= false;
00394         }
00395     }
00396     
00397     if (split) {
00398         setStatus(SPLIT);
00399         _mGFpair->kplanelist.push_back(_mGFpair->_mGFsegment->getKPlane());
00400     } else {
00401         _mGFpair->step(klayer);
00402     }
00403     _mGFpair->setStatus(FOUND);
00404     
00405 }
00406 //#########################################################################
00407 void GFpair::stepSplit(int klayer)
00408 //#########################################################################
00409 {
00410     _mGFbest->step(klayer);
00411     _mGFpair->step(klayer);
00412     
00413     if (_mGFbest->status() != FOUND || _mGFpair->status() != FOUND) return;
00414     
00415  //   if (_mGFbest->numDataPoints()==0 || _mGFpair->numDataPoints() == 0) GFcontrol::error++;
00416     
00417     int indexhit1 = _mGFbest->lastKPlane().getIDHit();
00418     int indexhit2 = _mGFpair->lastKPlane().getIDHit();
00419     if ( indexhit1 != indexhit2) return;
00420     
00421     selfishStepSplit(klayer);
00422 }
00423 
00424 //#########################################################################
00425 void GFpair::selfishStepSplit(int klayer)
00426 //#########################################################################
00427 {  
00428 
00429 //    if (_mGFbest->numDataPoints() <2 || _mGFpair->numDataPoints()<2) GFcontrol::error++;
00430     
00431     removeWorseStep(_mGFbest,_mGFpair);
00432     GFtrack* _GFwiner;
00433     GFtrack* _GFloser;
00434     if (_mGFpair->status() == EMPTY) {
00435         _GFwiner = _mGFbest;
00436         _GFloser = _mGFpair;
00437     } else {
00438         _GFwiner = _mGFpair;
00439         _GFloser = _mGFbest;
00440     }
00441     
00442     // flag the winer hit - remove the loser - next tray - unflag again
00443     
00444     int indexhit = _GFwiner->lastKPlane().getIDHit();
00445 
00446     GFtutor::_DATA->flagHit(m_axis, indexhit);
00447 
00448     _GFloser->step(klayer);
00449     GFtutor::_DATA->unflagHit(m_axis, indexhit);
00450     
00451     if (_GFloser->status() == EMPTY && allowedShareHit(_GFwiner)) {
00452         _GFloser->step(klayer);
00453         double sigma = _GFloser->lastKPlane().getSigma(KalHit::PRED);
00454         if (sigma > GFcontrol::maxSigmasSharedHits*m_sigmaCut) {
00455             _GFloser->removeStep();
00456         } 
00457     } 
00458     
00459 }
00460 //#########################################################################
00461 void GFpair::removeWorseStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2) 
00462 //#########################################################################
00463 {
00464     if (_GFtrack1->status() != FOUND || _GFtrack2->status() != FOUND) GFcontrol::error++;
00465     
00466     //unused:   int indexhit1 = -1;
00467     //unused:   int indexhit2 = -1;
00468     
00469     int kplane1 = _GFtrack1->previousKPlane().getIDPlane();
00470     int kplane2 = _GFtrack2->previousKPlane().getIDPlane();
00471     
00472     // the closest track first
00473 //    if (kplane1 != kplane2) 
00474 //    {
00475 //          if (kplane1 > kplane2) _GFtrack2->removeStep();
00476 //          else  _GFtrack1->removeStep();
00477 //          return;
00478 //    }
00479     
00480     int nhits1 = _GFtrack1->numDataPoints();
00481     int nhits2 = _GFtrack2->numDataPoints();
00482     
00483     double chi1 = _GFtrack1->lastKPlane().getDeltaChiSq(KalHit::FIT);
00484     double chi2 = _GFtrack2->lastKPlane().getDeltaChiSq(KalHit::FIT);
00485     
00486     if (nhits1 < 2 || nhits2 < 2) {
00487         if (_GFtrack1->_mGFsegment->accept() || 
00488             _GFtrack2->_mGFsegment->accept()  ) {
00489             if (_GFtrack1->_mGFsegment->accept()) chi1 = _GFtrack1->_mGFsegment->chiGFSq();
00490             else chi1 = 1e6;
00491             if (_GFtrack2->_mGFsegment->accept()) chi2 = _GFtrack2->_mGFsegment->chiGFSq();
00492             else chi2 = 1e6;
00493         }
00494     } 
00495 
00496     if (nhits2 < nhits1 || chi2 >= chi1) {
00497         _GFtrack2->removeStep();
00498     } else {
00499         _GFtrack1->removeStep();
00500     }
00501 }
00502 //#########################################################################
00503 bool GFpair::allowedShareHit(const GFtrack* _GFtrack) const
00504 //#########################################################################
00505 {
00506     bool allow = false;
00507     
00508     if (_GFtrack->numDataPoints() == 0) GFcontrol::error++;
00509 
00510     // the size of the cluster is big enough for this slope
00511     int idhit = _GFtrack->lastKPlane().getIDHit();
00512     double slope = _GFtrack->lastKPlane().getHit(KalHit::PRED).getPar().getSlope();
00513     int value = 2; // more strips than the expected
00514     if (GFtutor::okClusterSize(m_axis,idhit,slope) >= value) allow = true;
00515     
00516     return allow;
00517 }
00518 //#########################################################################
00519 Vector GFpair::doDirection(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2, 
00520                            double& weight1, double& errorSlope)
00521 //#########################################################################
00522 {
00523     Vector dir(0.,0.,0.);
00524     weight1 = 0.;
00525     errorSlope = 0.;
00526     if (_GFtrk1->empty()) return dir;
00527     
00528     double xene = _GFtrk1->m_iniEnergy/(_GFtrk1->m_iniEnergy+_GFtrk2->m_iniEnergy);
00529     double slopeBest   = _GFtrk1->slope();
00530     double slopePair   = _GFtrk2->slope();
00531     double errorSQbest = _GFtrk1->errorSlopeAtVertex();
00532     errorSQbest *= errorSQbest;
00533     double errorSQpair = _GFtrk2->errorSlopeAtVertex();
00534     errorSQpair *= errorSQpair;
00535     double errorSQ     = errorSQbest*errorSQpair/((1.-xene)*errorSQbest+xene*errorSQpair);
00536     
00537     weight1     = xene*errorSQ/errorSQbest;
00538     double slope = errorSQ*(xene*(slopeBest/errorSQbest)+(1.-xene)*(slopePair/errorSQpair));
00539     errorSlope  = sqrt(errorSQ);
00540     
00541     double factor = -1;
00542     double slopex = 0;
00543     double slopey =0;
00544     if (_GFtrk1->getAxis() == SiCluster::X) slopex = slope;
00545     else slopey = slope;
00546     
00547     dir = Vector(factor*slopex,factor*slopey,factor).unit();
00548     return dir;
00549 } 
00550 //########################################################
00551 Vector GFpair::doDirection(double& xFactor)
00552 //########################################################
00553 {
00554     Vector dir(0.,0.,0.);
00555     float wt1x = 1./(_mGFbest->chiSquare() + .01);
00556     float wt2x = 1./(3.*(_mGFpair->chiSquare() + .01));
00557     if (wt2x > wt1x) wt2x = wt1x;
00558     xFactor = wt1x/(wt1x + wt2x);
00559     Vector t1 = _mGFbest->direction();
00560     Vector t2 = _mGFpair->direction();
00561     //unused: float sinDLT12 = sqrt(std::max(0.0,(1. - sqr(t1*t2))));
00562     // if(sinDLT12 < .010/m_iniEnergy) { // only trust small opening angle paris
00563     double aveSlope = xFactor*_mGFbest->slope() + (1.-xFactor)*_mGFpair->slope();
00564     double slopeX = 0.;
00565     double slopeY = 0.;
00566     if (m_axis == SiCluster::X) slopeX = aveSlope;
00567     else slopeY = aveSlope;
00568     double factor = -1.;
00569     dir = Vector(factor*slopeX,factor*slopeY,factor).unit();
00570     // }  else dir = _mGFbest->direction();
00571     return dir;
00572 }
00573 //########################################################
00574 Vector GFpair::doDirectionXene(double xene, double& weight)
00575 //########################################################
00576 {
00577     Vector dir(0.,0.,0.);
00578     double x3 = xene*xene*xene;
00579     double ix3 = (1.-xene)*(1.-xene)*(1.-xene);
00580     float wt1x = x3/(x3+ix3);
00581     float wt2x = ix3/(x3+ix3);
00582     Vector t1 = _mGFbest->direction();
00583     Vector t2 = _mGFpair->direction();
00584     double aveSlope = wt1x*_mGFbest->slope() + wt2x*_mGFpair->slope();
00585     double slopeX = 0.;
00586     double slopeY = 0.;
00587     if (m_axis == SiCluster::X) slopeX = aveSlope;
00588     else slopeY = aveSlope;
00589     double factor = -1.;
00590     dir = Vector(factor*slopeX,factor*slopeY,factor).unit();
00591 
00592     weight = wt1x;
00593     return dir;
00594 }
00595 //#########################################################################
00596 double GFpair::doEnergy(const GFtrack* _GFtrk1, const GFtrack* _GFtrk2)
00597 //#########################################################################
00598 {
00599     double ene1 = _GFtrk1->RCenergy();
00600         if (m_iniEnergy < ene1) ene1 = m_iniEnergy;
00601     double ene2 = m_iniEnergy - ene1;
00602     if (!_GFtrk2->empty()) {
00603                 ene2 = _GFtrk2->RCenergy();
00604                 if (m_iniEnergy < ene2) ene2 = m_iniEnergy;
00605         }
00606     double x1 = ene1/m_iniEnergy;
00607     double x2 = 1.-ene2/m_iniEnergy;
00608     double x = 0.5*(x1+x2);
00609     return x;
00610 }
00611 
00612 //#########################################################################
00613 void GFpair::decideBest()
00614 //#########################################################################
00615 {
00616     // already decided!
00617     if (m_decideBest) return;
00618     
00619     double qbest = _mGFbest->doQbest();
00620     double qpair = _mGFpair->doQbest();
00621     
00622     if (qpair > qbest) swap();
00623     
00624     m_decideBest = true;
00625 }
00626 //#########################################################################
00627 void GFpair::swap()
00628 //#########################################################################
00629 {
00630     
00631     double eneBest = _mGFbest->m_iniEnergy;
00632     double enePair = _mGFpair->m_iniEnergy;
00633     double cutBest = _mGFbest->m_sigmaCut;
00634     double cutPair = _mGFbest->m_sigmaCut;
00635     
00636     _mGFalive = _mGFbest;
00637     _mGFbest  = _mGFpair;
00638     _mGFpair  = _mGFalive;
00639     
00640     _mGFbest->setIniEnergy(eneBest);
00641     _mGFpair->setIniEnergy(enePair);
00642     _mGFbest->m_sigmaCut = cutPair;
00643     _mGFpair->m_sigmaCut = cutBest;
00644     
00645     _mGFalive = (_mGFbest->m_alive? _mGFbest:_mGFpair);
00646     
00647 }
00648 
00649 //#########################################################################
00650 void GFpair::setIniEnergy()
00651 //#########################################################################
00652 {
00653     // set The real energies before the Fit!!
00654     _mGFbest->setIniEnergy(m_xEne*m_iniEnergy);
00655     _mGFpair->setIniEnergy((1.-m_xEne)*m_iniEnergy);
00656 }
00657 //-------------------------------------------------------------------------
00658 //
00659 //   Gamma
00660 //
00661 //-------------------------------------------------------------------------
00662 //#########################################################################
00663 GFgamma::GFgamma(double xene, 
00664                  double sigmaCut,
00665                  double energy, 
00666                  int ist, 
00667                  const Ray& testRay) : GFbase(sigmaCut,energy,ist,testRay),
00668                  m_xEne(xene)
00669                  //##########################################################################
00670 {
00671     
00672     GFdata::ini();
00673     
00674     m_connect = GFtutor::CONTROL_connectGFpair;
00675     m_associate = false;
00676     m_patternSwap = false;
00677     construct();
00678     if (!m_conflictPattern) return;
00679     
00680     clear();
00681     delete _mXpair;
00682     delete _mYpair;
00683     _mXpair = 0;
00684     _mYpair = 0;
00685     m_associate = true;
00686     // m_patternSwap;
00687     construct();
00688 }
00689 //#########################################################################
00690 void GFgamma::flagAllHits(int iflag)
00691 //#########################################################################
00692 {
00693     _mXpair->flagAllHits(iflag);
00694     _mYpair->flagAllHits(iflag);
00695 }
00696 //#########################################################################
00697 void GFgamma::unFlagAllHits()
00698 //#########################################################################
00699 {
00700     _mXpair->unFlagAllHits();
00701     _mYpair->unFlagAllHits();
00702 }
00703 //#########################################################################
00704 bool GFgamma::empty() const
00705 //##########################################################################
00706 {
00707     bool empty = GFdata::empty();
00708     empty = empty || _mXpair->empty() || _mYpair->empty();
00709     return empty;
00710 }
00711 
00712 //#########################################################################
00713 bool GFgamma::accept(const GFdata& pData1, const GFdata& pData2)
00714 //##########################################################################
00715 {
00716     bool ok = false; 
00717     
00718     int ilayerXfirst = pData1.firstLayer();
00719     int ilayerYfirst = pData2.firstLayer();
00720     if (abs(ilayerXfirst-ilayerYfirst) > GFcontrol::maxConsecutiveGaps) return ok;
00721     
00722     int iXtower = pData1.tower();
00723     int iYtower = pData2.tower();
00724     if (ilayerXfirst == ilayerYfirst) {
00725         if (iXtower == iYtower) ok = true;
00726         else ok = false;
00727     }
00728     
00729     return ok;
00730 }
00731 //#########################################################################
00732 void GFgamma::clear()
00733 //##########################################################################
00734 {
00735     m_status = TOGETHER;
00736     m_together = 0;
00737     m_split = 0;
00738     m_one = 0;
00739     
00740     _mXpair->clear();
00741     _mYpair->clear();
00742     
00743     GFdata::ini();
00744     setAlive();
00745     
00746 }
00747 
00748 
00749 //#########################################################################
00750 double GFgamma::Qbest()
00751 //##########################################################################
00752 {
00753     if (empty()) return m_quality;
00754     return _mXpair->_mGFbest->Qbest()+_mYpair->_mGFbest->Qbest();
00755 }
00756 //#########################################################################
00757 bool GFgamma::accept() const
00758 //##########################################################################
00759 {
00760     bool accept = false;
00761     if (empty()) return accept;
00762     
00763     accept = true;
00764     if (GFtutor::CUT_veto) accept = !veto();
00765     
00766     return accept;
00767 }
00768 
00769 //########################################################
00770 void GFgamma::writeOut(MsgStream& log) const
00771 //########################################################
00772 {
00773     // kludge to avoid warning messages from egcs
00774     int status = m_status;
00775     log << MSG::DEBUG << " --- GFgamma::writeOut --- " << endreq;
00776     log << MSG::DEBUG << " planes together " << numTogether() <<endreq;
00777     log << MSG::DEBUG << " planes split    " << numSplit() <<endreq;
00778     log << MSG::DEBUG << " planes one      " << numOne() << endreq;
00779     log << MSG::DEBUG << " last Status     " << status << endreq; 
00780     
00781 //    GFdata::writeOut(log);
00782     log << MSG::DEBUG << " --> Xpair : " <<endreq;
00783     _mXpair->writeOut(log);
00784     log << MSG::DEBUG << " --> Ypair : " <<endreq;
00785     _mYpair->writeOut(log);
00786     
00787 }
00788 
00789 //########################################################
00790 void GFgamma::draw(gui::DisplayRep& v) 
00791 //########################################################
00792 {
00793         _mXpair->draw(v);
00794         _mYpair->draw(v);
00795 
00796     // draw reconstructed gamma
00797     v.setColor("yellow");
00798     v.moveTo(this->vertex());
00799     v.lineTo(this->vertex()-300.*this->direction());
00800     v.setColor("black");
00801 
00802 
00803 }
00804 //#########################################################################
00805 bool GFgamma::veto() const
00806 //##########################################################################
00807 {
00808     bool veto = false;
00809     int ixhit;
00810     int iyhit;
00811     
00812     double sigma;
00813     // Veto for both views
00814     veto = _mXpair->_mGFbest->veto(ixhit,sigma) && _mYpair->_mGFbest->veto(iyhit,sigma);
00815     if (veto && (GFtutor::_DATA->getHit(SiCluster::X,ixhit)->tower() != 
00816                 GFtutor::_DATA->getHit(SiCluster::Y,iyhit)->tower())) veto = false;
00817     
00818     return veto; 
00819 }
00820 
00821 //#########################################################################
00822 Point GFgamma::getFirstHit() const
00823 //##########################################################################
00824 {
00825     Point firstHit;
00826     double zx = _mXpair->_mGFbest->vertex().z();
00827     double zy = _mYpair->_mGFbest->vertex().z();
00828     double xx = _mXpair->_mGFbest->vertex().x();
00829     double yy = _mYpair->_mGFbest->vertex().y();
00830     
00831     if (zy > zx) xx = _mXpair->_mGFbest->position(zy-zx);
00832     else yy = _mYpair->_mGFbest->position(zx-zy);
00833     
00834     double zz = (zy > zx? zy : zx);
00835     
00836     firstHit = Point(xx,yy,zz);
00837     return firstHit;
00838 }
00839 
00840 //-------------------------------------------------------------------------
00841 //   Gamma - Private
00842 //-------------------------------------------------------------------------
00843 
00844 //#########################################################################
00845 void GFgamma::ini() 
00846 //#########################################################################
00847 {
00848     
00849     _mXpair = 0;
00850     _mYpair = 0;
00851     _mXpair = new GFpair(m_xEne, SiCluster::X, m_sigmaCut, 
00852         m_iniEnergy, m_iniLayer, Ray(m_inVertex,m_inDirection), false);
00853     _mYpair = new GFpair(m_xEne, SiCluster::Y, m_sigmaCut, 
00854         m_iniEnergy, m_iniLayer, Ray(m_inVertex,m_inDirection), false);
00855     
00856     // controls
00857     m_fixTopology = false;
00858     m_decideBest = false;
00859     m_patternSwap = false;
00860     m_conflictPattern = false;
00861     m_swapDone = false;
00862     
00863     // status
00864     setDecideBest(m_decideBest);
00865     m_status = TOGETHER;
00866     
00867     //contability
00868     m_together = 0;
00869     m_split = 0;
00870     m_one = 0;
00871     
00872     setAlive();
00873     GFdata::ini();
00874     
00875 }
00876 
00877 //#########################################################################
00878 void GFgamma::step(int kplane) 
00879 //#########################################################################
00880 {
00881     if (!m_alive) return;
00882     
00883     _mXpair->step(kplane);
00884     _mYpair->step(kplane);
00885     
00886     m_status = newStatus();
00887     if (m_connect) {
00888         connectStep();
00889         if (!m_fixTopology) topologyStep();
00890 
00891         if (m_associate) associateStep();
00892     }
00893 }
00894 
00895 //#########################################################################
00896 void GFgamma::anastep(int kplane) 
00897 //#########################################################################
00898 {
00899     if (!m_alive) return;
00900 
00901     _mXpair->anastep(kplane);
00902     _mYpair->anastep(kplane);
00903 
00904     contability(kplane);
00905     if (m_associate) associateAnaStep();
00906     if (end()) {
00907         kill();
00908     }
00909 }
00910 
00911 //#########################################################################
00912 void GFgamma::fit()
00913 //##########################################################################
00914 {
00915     GFdata::ini();
00916     
00917     //  if (m_patternError) return;
00918     if (!m_decideBest) decideBest();
00919     
00920     _mXpair->fit();
00921     _mYpair->fit();
00922     if (_mYpair->empty()||_mXpair->empty()) return;
00923     
00924     loadGFdata();
00925     associateFit();
00926     //  if (m_associate) associateFit();
00927 }
00928 
00929 //#########################################################################
00930 bool GFgamma::end() const
00931 //#########################################################################
00932 {
00933     bool end = !m_alive;
00934     if (!_mXpair->m_alive || !_mYpair->m_alive) end = true;
00935     return end;
00936 }
00937 
00938 //#########################################################################
00939 void GFgamma::kill() 
00940 //#########################################################################
00941 {
00942     m_alive = false;
00943     _mXpair->kill();
00944     _mYpair->kill();
00945 }
00946 
00947 //#########################################################################
00948 void GFgamma::setAlive() 
00949 //#########################################################################
00950 {
00951     m_alive = true;
00952     _mXpair->setAlive();
00953     _mYpair->setAlive();
00954 }
00955 
00956 //#########################################################################
00957 void GFgamma::contability(int kplane) 
00958 //#########################################################################
00959 {
00960     if (m_status == TOGETHER) m_together++;
00961     else if (m_status == SPLIT) m_split++;
00962     else if (m_status == ONE) m_one++; 
00963 }
00964 
00965 //#########################################################################
00966 void GFgamma::loadGFdata()
00967 //##########################################################################
00968 {
00969     m_quality = _mXpair->Q() + _mYpair->Q();
00970     
00971     m_direction = GFdata::doDirection(_mXpair->direction(),_mYpair->direction());
00972     
00973     Ray XRay = Ray(_mXpair->vertex(),_mXpair->direction());
00974     Ray YRay = Ray(_mYpair->vertex(),_mYpair->direction());
00975     m_vertex=GFbase::doVertex( XRay, YRay);
00976     
00977     m_RCenergy = 0.5*(_mXpair->RCenergy()+_mYpair->RCenergy());
00978     
00979     int ixlayer = _mXpair->firstLayer();
00980     int iylayer = _mYpair->firstLayer();
00981     m_firstLayer = (ixlayer < iylayer? ixlayer:iylayer);
00982     
00983     m_nhits = _mXpair->nhits()+_mYpair->nhits();
00984     
00985     m_itower = ( ixlayer < iylayer? _mXpair->tower() : _mYpair->tower());
00986 }
00987 //#########################################################################
00988 void GFgamma::construct()
00989 //##########################################################################
00990 {
00991     ini();
00992 
00993     doit();
00994 
00995     fit();
00996 
00997     if (empty()) clear();
00998 }
00999 
01000 //#########################################################################
01001 GFbase::StatusPair GFgamma::newStatus() 
01002 //#########################################################################
01003 {
01004     StatusPair Xstatus = _mXpair->status();
01005     StatusPair Ystatus = _mYpair->status();
01006     StatusPair status = Xstatus;
01007     
01008     if (Xstatus == Ystatus) status = Xstatus;
01009     else if (Xstatus == DONE || Ystatus == DONE) status = DONE; 
01010     else if (Xstatus == SPLIT || Ystatus == SPLIT) status = SPLIT;
01011     else if (Xstatus == TOGETHER || Ystatus == TOGETHER) status = TOGETHER;
01012     return status;
01013 }
01014 //#########################################################################
01015 void GFgamma::connectStep() 
01016 //#########################################################################
01017 {
01018     if (m_status == TOGETHER) {
01019 
01020         if (!GFparticle::sameTower(_mXpair->_mGFbest, _mYpair->_mGFbest)) {
01021             
01022             bool done = GFparticle::removeWorseStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
01023             if (!done) m_conflictPattern = true;
01024         }
01025     }
01026     
01027     if (!m_associate) {
01028         if (m_status == SPLIT) {
01029 
01030             bool bestbest = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFbest);
01031             bool pairpair = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFpair);
01032             if (!bestbest || !pairpair) {
01033 
01034                 bool bestpair = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFpair);
01035                 bool pairbest = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFbest);
01036                 if (!bestpair || !pairbest) m_conflictPattern = true;
01037             }
01038         }
01039         if (m_status == ONE) {
01040 
01041             if (!GFparticle::sameTower(_mXpair->_mGFalive,_mYpair->_mGFalive)) {
01042 
01043                 bool done = GFparticle::removeWorseStep(_mXpair->_mGFalive,_mYpair->_mGFalive);
01044                 if (!done) m_conflictPattern = true;
01045             }
01046         }
01047     }
01048     
01049 }
01050 
01051 //#########################################################################
01052 void GFgamma::topologyStep() 
01053 //#########################################################################
01054 {
01055     if (m_status == SPLIT) {
01056         bool bestbest = crossingTowers(_mXpair->_mGFbest,_mYpair->_mGFbest,
01057             _mXpair->_mGFpair,_mYpair->_mGFpair);
01058         bool bestpair = crossingTowers(_mXpair->_mGFbest,_mYpair->_mGFpair,
01059             _mXpair->_mGFpair,_mYpair->_mGFbest);
01060         if (bestbest && bestpair) GFcontrol::error++;
01061         else if (bestbest || bestpair) {
01062             m_fixTopology = true;
01063             if (!m_associate) {
01064                 m_associate = true;
01065                 m_patternSwap = bestpair;
01066             }
01067         }
01068     } else if (m_status == ONE) {
01069         if (_mXpair->_mGFalive->status() == FOUND && 
01070             _mYpair->_mGFalive->status() == FOUND) {
01071             m_fixTopology = true;
01072             if (!m_associate) {
01073                 m_associate = true;
01074                 if ((_mXpair->_mGFbest->alive() && _mYpair->_mGFbest->alive()) ||
01075                     (_mXpair->_mGFpair->alive() && _mYpair->_mGFpair->alive()) )
01076                     m_patternSwap = false;
01077                 else m_patternSwap = true;
01078             }
01079         }
01080     }
01081 }
01082 //#########################################################################
01083 void GFgamma::associateStep() 
01084 //#########################################################################
01085 {
01086     associateStatus(m_status);
01087     
01088     if (m_status == SPLIT) {
01089         // TOGETHER already tested in connectStep;
01090         bool bestbest = GFparticle::sameTower(_mXpair->_mGFbest,_mYpair->_mGFbest);
01091         bool pairpair = GFparticle::sameTower(_mXpair->_mGFpair,_mYpair->_mGFpair);
01092         if (!bestbest || !pairpair) {
01093             bool done = false;
01094             // spetial case - first plane
01095 
01096             if (!bestbest) done = GFparticle::removeWorseStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
01097             if (!pairpair) done = GFparticle::removeWorseStep(_mXpair->_mGFpair,_mYpair->_mGFpair);
01098             if (!done) m_conflictPattern = true;
01099         }
01100     } else if (m_status == ONE) {
01101         if (!GFparticle::sameTower(_mXpair->_mGFalive, _mYpair->_mGFalive)) { 
01102 
01103             bool done = GFparticle::removeWorseStep(_mXpair->_mGFalive,_mYpair->_mGFalive);
01104             if (!done) m_conflictPattern = true;
01105         }
01106     }
01107 }
01108 
01109 //#########################################################################
01110 void GFgamma::associateStatus(StatusPair status) 
01111 //#########################################################################
01112 {       
01113     StatusPair Xstatus = _mXpair->status();
01114     StatusPair Ystatus = _mYpair->status();
01115     if (m_patternSwap && !m_swapDone && m_status != TOGETHER) {
01116         //      if (!_mXpair->_mGFpair->m_alive && !_mYpair->_mGFpair->m_alive) GFcontrol::error++;
01117         if (Xstatus == SPLIT || _mXpair->numSplit()>0) _mXpair->swap();
01118         else if (Ystatus == SPLIT || _mYpair->numSplit()>0) _mYpair->swap();
01119         else GFcontrol::error++;
01120         m_swapDone = true; // reset the pattern - now again best-best.
01121     }
01122     
01123     _mXpair->setStatus(status);
01124     _mYpair->setStatus(status);
01125     
01126     m_status = status;
01127 }
01128 //#########################################################################
01129 bool GFgamma::crossingTowers(const GFtrack* _GFtrkX1, const GFtrack* _GFtrkY1,
01130                              const GFtrack* _GFtrkX2, const GFtrack* _GFtrkY2) 
01131                              //#########################################################################
01132 {
01133     // The two tracks are in different towers
01134     bool fix = false;
01135     
01136     int nX1hits = _GFtrkX1->numDataPoints();
01137     int nY1hits = _GFtrkY1->numDataPoints();
01138     if ( nX1hits < 1 || nY1hits < 1 ) return fix;
01139     
01140     int nX2hits = _GFtrkX2->numDataPoints();
01141     int nY2hits = _GFtrkY2->numDataPoints();
01142     if ( nX2hits < 1 || nY2hits < 1 ) return fix;
01143     
01144     if (_GFtrkX1->status() != FOUND || _GFtrkY1->status() != FOUND ||
01145         _GFtrkX2->status() != FOUND || _GFtrkY2->status() != FOUND ) return fix;
01146     
01147     int iX1tower = _GFtrkX1->kplanelist[nX1hits-1].getIDTower();
01148     int iY1tower = _GFtrkY1->kplanelist[nY1hits-1].getIDTower();
01149     int iX2tower = _GFtrkX2->kplanelist[nX2hits-1].getIDTower();
01150     int iY2tower = _GFtrkY2->kplanelist[nY2hits-1].getIDTower();
01151     
01152     if ((iX1tower == iY1tower) &&
01153         (iX2tower == iY2tower) &&
01154         (iX1tower != iX2tower)) fix = true;
01155     
01156     return fix;
01157 }
01158 
01159 
01160 //#########################################################################
01161 void GFgamma::associateAnaStep() 
01162 //#########################################################################
01163 {
01164     associateAnaStep(_mXpair->_mGFbest,_mYpair->_mGFbest);
01165     associateAnaStep(_mYpair->_mGFbest,_mXpair->_mGFbest);
01166     associateAnaStep(_mXpair->_mGFpair,_mYpair->_mGFpair);
01167     associateAnaStep(_mYpair->_mGFpair,_mXpair->_mGFpair);
01168 }
01169 
01170 //#########################################################################
01171 void GFgamma::associateAnaStep(GFtrack* _GFtrack1, GFtrack* _GFtrack2) 
01172 //#########################################################################
01173 {
01174     if (_GFtrack1->numDataPoints() == 0) return;
01175     if (_GFtrack1->status() != FOUND) return;
01176     
01177     if (!_GFtrack2->m_alive) {
01178                 if (_GFtrack2->numDataPoints() < GFcontrol::minSegmentHits) _GFtrack1->clear();
01179         _GFtrack1->kill();
01180     } else {
01181         _GFtrack1->associateOrthStep(_GFtrack2);
01182     }
01183 }
01184 //#########################################################################
01185 void GFgamma::associateFit() 
01186 //#########################################################################
01187 {
01188     bool fix = m_fixTopology;
01189     KalHit::TYPE typ = KalHit::SMOOTH;
01190     _mXpair->_mGFbest->associateOrthGFtrack(_mYpair->_mGFbest, fix, typ);
01191     _mYpair->_mGFbest->associateOrthGFtrack(_mXpair->_mGFbest, fix, typ);
01192     _mXpair->_mGFpair->associateOrthGFtrack(_mYpair->_mGFpair, fix, typ);
01193     _mYpair->_mGFpair->associateOrthGFtrack(_mXpair->_mGFpair, fix, typ);
01194 }
01195 
01196 //#########################################################################
01197 void GFgamma::setDecideBest(bool decideBest)
01198 //##########################################################################
01199 {
01200     m_decideBest = decideBest;
01201     _mXpair->setDecideBest(m_decideBest);
01202     _mYpair->setDecideBest(m_decideBest);
01203 }
01204 //#########################################################################
01205 void GFgamma::decideBest()
01206 //##########################################################################
01207 {
01208     
01209     if (!m_associate) return;
01210     
01211     setDecideBest(true);
01212     
01213     double qbest = _mXpair->_mGFbest->doQbest();
01214     qbest += _mYpair->_mGFbest->doQbest();
01215     
01216     double qpair = _mXpair->_mGFpair->doQbest();
01217     qpair += _mYpair->_mGFpair->doQbest();
01218     
01219     if (qbest < qpair) {
01220         _mXpair->swap();
01221         _mYpair->swap();
01222     }
01223     
01224 }