GFcandidates.cpp Source File

00001 #include "TkrRecon/GFcandidates.h"
00002 
00003 //###########################################################
00004 GFcandidates::GFcandidates(enum GFcandidates::type t, double ene, double sigmaCut,
00005                                                    Point Pend, Point Pini):m_type(t),
00006                                                    m_eneCandidate(ene),m_sigmaCut(sigmaCut),m_Pend(Pend),m_Pini(Pini)
00007 //###########################################################
00008 {
00009         ini();
00010 
00011         bool okX = findSeedCandidates(m_Xcandidates,m_seedtype,SiCluster::X);
00012         bool okY = findSeedCandidates(m_Ycandidates,m_seedtype,SiCluster::Y);
00013 
00014         if (m_type != m_seedtype) findCandidates();
00015 
00016 }
00017 
00018 //###########################################################
00019 void GFcandidates::clear()
00020 //###########################################################
00021 {
00022         m_Xcandidates.clear();
00023         m_Ycandidates.clear();
00024         m_candidates.clear();
00025 }
00026 
00027 //###########################################################
00028 GFdata GFcandidates::GFconstructor(enum GFcandidates::type type, double ene, double sigmaCut,
00029                                                                    int ilayer,const Ray testRay, SiCluster::view axis)
00030 //###########################################################
00031 {
00032     GFdata data;
00033 
00034     if (type == GFcandidates::PARTICLE) {
00035         GFparticle* _par = new GFparticle(sigmaCut,     ene, ilayer, testRay);
00036         if (!_par->empty() && _par->accept()) data = _par->getGFdata();
00037         delete _par;
00038     } else if (type == GFcandidates::GAMMA) {
00039         GFgamma* _gamma = new GFgamma(GFcontrol::FEne, sigmaCut, ene, ilayer, testRay);
00040         if (!_gamma->empty() && _gamma->accept()) {
00041             data = _gamma->getGFdata();
00042         }
00043         delete _gamma;
00044     } else if (type == GFcandidates::TRACK) {
00045         GFtrack* _track = new GFtrack(axis, sigmaCut, ene, ilayer, testRay);
00046         if (!_track->empty() && _track->accept()) data = _track->getGFdata();
00047         delete _track;
00048     } else if (type == GFcandidates::PAIR) {
00049         GFpair* _pair = new GFpair(GFcontrol::FEne, axis, sigmaCut, ene, ilayer, testRay);
00050         if (!_pair->empty()) {
00051             if (_pair->accept()) {
00052                 data = _pair->getBest()->getGFdata();
00053             }
00054         }
00055         delete _pair;
00056     }
00057 
00058     return data;
00059 }
00060 //-----------  Private drivers  ----------------------------- 
00061 
00062 //###########################################################
00063 void GFcandidates::ini()
00064 //###########################################################
00065 {
00066         clear();
00067         if (m_type == GAMMA) m_seedtype = PAIR;
00068         if (m_type == PARTICLE) m_seedtype = TRACK;
00069 }
00070 
00071 //###########################################################
00072 bool  GFcandidates::findCandidates(std::vector<GFdata>& candidates,
00073                                                                    const GFdata& Xcandidate, 
00074                                                                    const GFdata& Ycandidate,
00075                                                                    double ene,
00076                                                                    enum GFcandidates::type typ)
00077 //###########################################################
00078 {
00079     bool ok = false;
00080     int naccepted = 0;
00081 
00082     int iniLayer = (Xcandidate.firstLayer() < Ycandidate.firstLayer()?
00083         Xcandidate.firstLayer() : Ycandidate.firstLayer());
00084     int lastLayer = (Xcandidate.firstLayer()< Ycandidate.firstLayer()?
00085         Ycandidate.firstLayer() : Xcandidate.firstLayer());
00086 
00087     if (lastLayer - iniLayer >= GFcontrol::maxConsecutiveGaps) return ok;
00088 
00089     Point ver =GFdata::doVertex(Xcandidate.ray(),Ycandidate.ray());
00090     Vector dir = GFdata::doDirection(Xcandidate.direction(),Ycandidate.direction());
00091     Ray testRay(ver, dir);
00092 
00093     for (int ilayer = iniLayer ; ilayer <= lastLayer; ilayer++) {
00094 
00095         GFdata candidateGFdata = GFconstructor(typ, m_eneCandidate, m_sigmaCut, ilayer, testRay);
00096         if (candidateGFdata.Q() > GFcontrol::minQ) {
00097             ok = true;
00098             naccepted++;
00099             incorporate(candidates, candidateGFdata);
00100         }
00101     }
00102 
00103     return ok;
00104 }
00105 //------------- Utilities -----------------------------------
00106 
00107 //###########################################################
00108 bool GFcandidates::findSeedCandidates(std::vector<GFdata>& candidates, 
00109                                                                           GFcandidates::type typ, SiCluster::view axis)
00110 //###########################################################
00111 {
00112     bool OK = false;
00113     for (int iplane = 0 ; iplane < GFtutor::numPlanes() - 2; iplane++){
00114         bool ok = findSeedCandidates(candidates, typ, axis, iplane);
00115         OK = OK || ok;
00116     }
00117     return OK;
00118 }
00119 
00120 //###########################################################
00121 bool GFcandidates::findSeedCandidates(std::vector<GFdata>& candidates, 
00122                                                                           GFcandidates::type typ, SiCluster::view axis,
00123                                                                           int ilayer, int itower)
00124 //###########################################################
00125 {
00126     //unused:   int nconstructed = 0;
00127     int naccepted = 0;
00128     bool ok = false;
00129 
00130         Point Pini = Point(0.,0.,0.);
00131         Point Pend = Point(0.,0.,0.);
00132 
00133         bool loopHits = true;
00134         if (m_Pini.mag() > 0.001) loopHits = false;
00135         
00136 
00137         std::vector<SiCluster*> hitList;
00138     int nhits = GFtutor::_DATA->nHits(axis,ilayer);
00139         if (nhits > 0) hitList = GFtutor::_DATA->getHits(axis,ilayer);
00140 
00141         bool end = false;
00142         if (nhits == 0) end = true;
00143 
00144         int ihit = 0;
00145         while (!end) {
00146         
00147                 if (loopHits) Pini = hitList[ihit]->position();
00148                 else Pini = m_Pini;
00149 
00150         Pend = createPend(axis, ilayer, Pini);
00151 
00152         Vector VDir(Pend.x()-Pini.x(),Pend.y()-Pini.y(),Pend.z()-Pini.z());
00153         Ray testRay = Ray(Pini, VDir.unit());
00154 
00155         GFdata candidateGFdata = GFconstructor(typ, m_eneCandidate, m_sigmaCut, ilayer, testRay, axis);
00156 
00157         if (candidateGFdata.Q() > GFcontrol::minQ) {
00158             ok = true;
00159             naccepted++;
00160             incorporate(candidates, candidateGFdata);
00161         }
00162 
00163                 if (!loopHits) end = true;
00164                 else ihit++;
00165                 if (ihit >= nhits) end = true;
00166     }
00167 
00168     if (naccepted > 0) ok = true;
00169     return ok;
00170 
00171 }
00172 //###########################################################
00173 void GFcandidates::incorporate(std::vector<GFdata>& pDatalist, const GFdata pData)
00174 //###########################################################
00175 {
00176     bool ienter = false;
00177     for (unsigned int i =0; i < pDatalist.size(); i++) {
00178         if (pData.Q() >= pDatalist[i].Q()) {
00179             pDatalist.insert(&pDatalist[i],pData);
00180             ienter = true;
00181         }
00182         if (ienter) break;
00183     }
00184     if (!ienter) pDatalist.push_back(pData);
00185     if (pDatalist.size()>GFcontrol::maxCandidates) pDatalist.pop_back();
00186 
00187 }
00188 
00189 //########################################################
00190 Point GFcandidates::createPend(SiCluster::view axis,int ilayer, const Point& Pini)
00191 //########################################################
00192 {
00193     double weight = GFcontrol::minEnergy/(3.*m_eneCandidate); // weight;
00194 
00195     Point PCal = m_Pend;
00196 
00197     if (m_eneCandidate < GFcontrol::minEnergy) weight = 1.;
00198     if (PCal.mag() == 0) weight = 1.;
00199         //double side = GFtutor::trayWidth();
00200     //I think the above is too big and causing problems... 
00201     //Look in a region below the current hit which is within a cone slightly 
00202     //larger than 45 degrees (arbitrary!) of the current hit.
00203     //double side = 2.5 * GFtutor::trayGap();
00204     //double side = 3.5 * GFtutor::trayGap();
00205     double side = 5.0 * GFtutor::trayGap();
00206     Point PTrk = GFtutor::_DATA -> meanHitInside(axis, ilayer+1,0.5*side, Pini);
00207 
00208     if (PTrk.mag() == 0.) 
00209     {
00210         side = 2 * side;
00211         PTrk = GFtutor::_DATA -> meanHitInside(axis, ilayer+2,0.5*side, Pini);
00212     }
00213         
00214     if (PTrk.mag() == 0.) weight = 0.;
00215     double x = (1.-weight)*PCal.x()+weight*PTrk.x();
00216     double y = (1.-weight)*PCal.y()+weight*PTrk.y();
00217     double z = (1.-weight)*PCal.z()+weight*PTrk.z();
00218 
00219     Point PRef(x,y,z);
00220 
00221     if (PRef.mag() != 0.) 
00222     {
00223         //Add an offset to the mean position to prevent deadlock when only two hits used
00224         if (axis == SiCluster::X ) 
00225         {
00226             //x += 0.5 * GFtutor::siResolution();
00227             y  = Pini.y();
00228         }
00229         else 
00230         {
00231             x  = Pini.x();
00232             //y += 0.5 * GFtutor::siResolution();
00233         }
00234     } else 
00235     {
00236         x = Pini.x();
00237         y = Pini.y();
00238         z = Pini.z()-GFtutor::trayGap();
00239     }
00240 
00241     PRef = Point(x,y,z);
00242 
00243     return PRef;
00244 }
00245 
00246 //###########################################################
00247 bool GFcandidates::findCandidates()
00248 //###########################################################
00249 {
00250     bool OK = false;
00251 
00252     m_candidates.clear();
00253     int ix = 0;
00254     for (; ix < m_Xcandidates.size(); ix++) {
00255         for (int iy = 0; iy < m_Ycandidates.size(); iy++) {
00256             GFdata Xcandidate = m_Xcandidates[ix];
00257             GFdata Ycandidate = m_Ycandidates[iy];
00258             bool ok = false;
00259             ok = findCandidates(m_candidates, 
00260                                 Xcandidate, Ycandidate, m_eneCandidate, m_type);
00261             OK = OK || ok;
00262         } // Y candidates
00263     } // X candidates
00264 
00265     if (OK)  return OK;
00266 
00267     // Force a PairFit with no veto
00268     bool save_veto = GFtutor::CUT_veto;
00269     std::vector<GFdata> candidates;
00270     candidates.clear();
00271     for (ix = 0 ; ix < m_Xcandidates.size(); ix++) {
00272         GFdata Xcandidate = m_Xcandidates[ix];
00273         candidates.clear();
00274         GFtutor::CUT_veto = false;
00275 
00276         findSeedCandidates(candidates, m_seedtype, SiCluster::Y, Xcandidate.firstLayer(), Xcandidate.tower());
00277         for (int iy =0 ; iy < candidates.size(); iy++) {
00278             GFtutor::CUT_veto = save_veto;
00279             GFdata Ycandidate = candidates[iy];
00280             bool ok = findCandidates(m_candidates, 
00281                                 Xcandidate,Ycandidate, m_eneCandidate, m_type);
00282             OK = OK || ok;
00283         } // Y candidates
00284     } // X candidates
00285 
00286     candidates.clear();
00287     for (int iy = 0 ; iy < m_Ycandidates.size(); iy++) {
00288         GFdata Ycandidate = m_Ycandidates[iy];
00289         candidates.clear();
00290         GFtutor::CUT_veto = false;
00291 
00292         findSeedCandidates(candidates, m_seedtype, SiCluster::X, Ycandidate.firstLayer(), Ycandidate.tower());
00293         for (ix =0 ; ix < candidates.size(); ix++) {
00294             GFtutor::CUT_veto = save_veto;
00295             GFdata Xcandidate = candidates[ix];
00296             bool ok = findCandidates(m_candidates, 
00297                                 Xcandidate,Ycandidate, m_eneCandidate, m_type);
00298             OK = OK || ok;
00299         } // Y candidates
00300     } // X candidates
00301     GFtutor::CUT_veto = save_veto;
00302     return OK;
00303 }