#include <FindTrackSeg.h>
Inheritance diagram for rpr::FindTrackSeg:
Public Member Functions | |
| FindTrackSeg (const char *version) | |
| void | Update (const cfg::Config &c) |
| virtual | ~FindTrackSeg () |
| jobc::Result | Reco (edm::EventHandle &evt) |
| void | DoView (std::vector< const recobase::PlaneCluster * > planeclust, std::vector< recobase::TrackBase > &track, geo::_plane_proj xychoice) |
Private Attributes | |
| short | fDetGeom |
| std::vector< float > | fGapThreshold |
| Allowed gap sizes in tracks. | |
| int | fMaxTracks |
| Max. tracks allowed per view. | |
| unsigned int | fMinTrackHits |
| Min. hits allowed to form a track. | |
| float | fHoughThreshold |
| Threshold to apply to Hough map. | |
| geo::Geometry * | fGeo |
| double | fInterval |
| Distance between planes. | |
|
|
Definition at line 31 of file FindTrackSeg.cxx. References jobc::Module::SetCfgVersion().
00031 : 00032 jobc::Module("FindTrackSeg"), 00033 fDetGeom(1), fGeo(0), fInterval(13.27) 00034 { 00035 this->SetCfgVersion(version); 00036 } |
|
|
Definition at line 40 of file FindTrackSeg.cxx.
00041 {
00042 }
|
|
||||||||||||||||
|
Definition at line 102 of file FindTrackSeg.cxx. References rpr::Hough::AddPoint(), recobase::TrackBase::Clear(), recobase::PlaneCluster::dT(), recobase::PlaneCluster::dZ(), recobase::TrackBase::Erase(), fGapThreshold, fHoughThreshold, rpr::Line::fIcept, rpr::Hough::FindLines(), fInterval, fMaxTracks, fMinTrackHits, rpr::Line::fScore, rpr::Line::fSlope, recobase::RecoPoint::GetPC(), rpr::Hough::MakeTree(), recobase::TrackBase::NRecoPt(), recobase::PlaneCluster::Plane(), recobase::TrackBase::Put(), recobase::TrackBase::RecoPt(), recobase::RecoPoint::Set(), rpr::Hough::SetThresh(), recobase::TrackBase::SetXY(), rpr::Hough::TestHit(), recobase::PlaneCluster::Tpos(), rpr::Hough::UpdateTree(), recobase::RecoPoint::X(), recobase::RecoPoint::Y(), recobase::RecoPoint::Z(), and recobase::PlaneCluster::Zpos(). Referenced by Reco().
00105 {
00106
00107 if(clust.size()<2) return;
00108
00109 vector<int> infit;
00110 infit.resize(clust.size());
00111 for (unsigned int i=0; i<infit.size(); ++i) infit[i] = 0;
00112
00113 Hough hough;
00114
00115 /* Hough Transform
00116 1. Fill point to theta vs rho Hough map
00117 2. Pick up maximum point in the map
00118 3. Find associate hits with the track
00119 4. remove isolate hits
00120 5. subtract hits from Hough map
00121
00122 end conditions:
00123 in 2. if score is below threshold
00124 in 3. if no hit was found
00125
00126 exceptions:
00127 in 5. if only one plane is in the track, pick up next maximum
00128 */
00129
00130 int nfound=0;
00131
00132 hough.SetThresh(fHoughThreshold);
00133
00134 float Zave = 0.;
00135 float Tave = 0.;
00136
00137 vector<const PlaneCluster*>::iterator itr = clust.begin();
00138 while(itr!=clust.end()) {
00139 Zave += (*itr)->Zpos();
00140 Tave += (*itr)->Tpos();
00141 ++itr;
00142 }
00143 Zave /= clust.size();
00144 Tave /= clust.size();
00145
00146 // Make Hough map
00147 // 1. Fill point to theta vs rho Hough map
00148
00149 for(unsigned int ih1=0; ih1<clust.size()-1; ih1++) {
00150 for(unsigned int ih2=ih1+1; ih2<clust.size(); ih2++) {
00151 // should be in different planes
00152 if(clust[ih1]->Plane()!=clust[ih2]->Plane()) {
00153
00154 // Use number of clusters instead of charge to be uniform in detector
00155 hough.AddPoint(clust[ih1]->Zpos()-Zave, clust[ih1]->dZ(),
00156 clust[ih2]->Zpos()-Zave, clust[ih2]->dZ(),
00157 clust[ih1]->Tpos()-Tave, clust[ih1]->dT(),
00158 clust[ih2]->Tpos()-Tave, clust[ih2]->dT(),
00159 1., ih1);
00160 }
00161 }
00162 }
00163 hough.MakeTree(fHoughThreshold);
00164
00165 TrackBase trk;
00166 RecoPoint pt;
00167 TVector3 pos;
00168
00169 // Search track candidates up to 100
00170 while(nfound<fMaxTracks) {
00171
00172 trk.Clear();
00173
00174 Line ln;
00175 double theta, rho, score;
00176
00177 // 2. Pick up maximum point in the map
00178 hough.FindLines(&theta, &rho, &score);
00179 if (score<fHoughThreshold) break;
00180
00181 double a, b;
00182 b = -cos(theta)/sin(theta);
00183 a = rho/sin(theta)+Tave-b*Zave;
00184 ln.fScore = score;
00185 ln.fIcept = a;
00186 ln.fSlope = b;
00187
00188 //use this version if you only want to consider hits which made up
00189 //the hough bin 3.
00190 /*
00191 int hit_to_add=-1;
00192 hit_to_add=hough.getHit();
00193 while (hit_to_add>-1){
00194 if(infit[hit_to_add]==0)
00195 {
00196 infit[hit_to_add]=nfound+1;
00197
00198 double x = (xychoice==geo::kX) ?
00199 ln.fIcept + clust[hit_to_add]->Zpos()*ln.fSlope : 0.;
00200 double y = (xychoice==geo::kY) ?
00201 ln.fIcept + clust[hit_to_add]->Zpos()*ln.fSlope : 0.;
00202 double z = clust[hit_to_add]->Zpos();
00203 pos.SetXYZ(x,y,z);
00204
00205 pt.Set(pos, clust[hit_to_add], 1.0);
00206
00207 trk.Put(pt);
00208 }
00209 hit_to_add=hough.getHit();
00210 }
00211 */
00212
00213 // 3. Find associate hits with the track
00214 for(unsigned int i=0; i<clust.size(); i++){
00215 if (infit[i]==0 && hough.TestHit(clust[i]->Zpos() - Zave,
00216 clust[i]->Tpos() - Tave,
00217 theta, rho,
00218 // 2.0*M_PI/180.0, 0.03)==true) {
00219 clust[i]->dT()*2.)==true) {
00220 infit[i] = nfound + 1;
00221
00222 double x =
00223 (xychoice==geo::kX) ? ln.fIcept + clust[i]->Zpos()*ln.fSlope : 0.;
00224 double y =
00225 (xychoice==geo::kY) ? ln.fIcept + clust[i]->Zpos()*ln.fSlope : 0.;
00226 double z = clust[i]->Zpos();
00227 pos.SetXYZ(x,y,z);
00228
00229 pt.Set(pos, clust[i], 1.0);
00230
00231 trk.Put(pt);
00232 }
00233 }
00235
00236 // Finish search if no hit is accosiated with track
00237 if(trk.NRecoPt()<=1) break;
00238
00239 // Only one hit is allowed in one plane.
00240 unsigned int i = 0;
00241 while(1){
00242 if(trk.RecoPt(i).GetPC()->Plane()==trk.RecoPt(i+1).GetPC()->Plane()) {
00243
00244 float dist1 =
00245 fabs(cos(theta)*(trk.RecoPt(i).GetPC()->Zpos() - Zave) +
00246 sin(theta)*(trk.RecoPt(i).GetPC()->Tpos() - Tave) -
00247 rho);
00248 float dist2 =
00249 fabs(cos(theta)*(trk.RecoPt(i+1).GetPC()->Zpos() - Zave) +
00250 sin(theta)*(trk.RecoPt(i+1).GetPC()->Tpos() - Tave) -
00251 rho);
00252
00253 if( dist1/trk.RecoPt(i).GetPC()->dT() >
00254 dist2/trk.RecoPt(i+1).GetPC()->dT() ) {
00255 for(unsigned int j=0; j<clust.size(); j++) {
00256 if(clust[j] == trk.RecoPt(i).GetPC()) infit[j] = 0;
00257 }
00258 trk.Erase(i);
00259 }
00260 else {
00261 for(unsigned int j=0; j<clust.size(); j++)
00262 if(clust[j] == trk.RecoPt(i+1).GetPC()) infit[j] = 0;
00263 trk.Erase(i+1);
00264 }
00265 }
00266 else {
00267 ++i;
00268 }
00269 if(i >= trk.NRecoPt()-1) break;
00270 }
00271
00272 // 4. Remove isolated hits
00273 float exgap = 1.0E10;
00274 unsigned int exhit = 0;
00275 while(1){
00276 if(trk.NRecoPt()<=1) break;
00277 // count number of hits before and after maximum gap
00278 int hitsaftergap = 0;
00279 unsigned int igap=0;
00280 float gap, maxgap = -1.;
00281
00282 for(unsigned int i=1; i<trk.NRecoPt(); i++) {
00283 gap = sqrt(pow(trk.RecoPt(i).X() - trk.RecoPt(i-1).X(),2)
00284 + pow(trk.RecoPt(i).Y() - trk.RecoPt(i-1).Y(),2)
00285 + pow(trk.RecoPt(i).Z() - trk.RecoPt(i-1).Z(),2));
00286 if(gap>maxgap && (gap<exgap || (gap==exgap && i>exhit))) {
00287 hitsaftergap=0;
00288 igap = i;
00289 maxgap = gap;
00290 }
00291 ++hitsaftergap;
00292 }
00293 exgap = maxgap;
00294 exhit = igap;
00295 // Gap condition: while loop until gap becomes smaller than the condition
00296 if(maxgap<fInterval*fGapThreshold[0]) break; // (# planes in X(Y) view only)
00297
00298 // Criteria to remove isolated hits dependes on the number of
00299 // hits after gap
00300 float threshold;
00301 switch(min(hitsaftergap, (int)trk.NRecoPt() - hitsaftergap)){
00302 case 1:
00303 threshold = fGapThreshold[0];
00304 break;
00305 case 2:
00306 threshold = fGapThreshold[1];
00307 break;
00308 default:
00309 threshold = fGapThreshold[2];
00310 break;
00311 }
00312
00313 if(maxgap>fInterval*threshold) {
00314 // Remove hits in shorter segment
00315 if((float)hitsaftergap<0.5*(float)trk.NRecoPt()) {
00316 for (unsigned int i=igap; i<trk.NRecoPt(); ++i) {
00317 for(unsigned int j=0; j<clust.size(); j++)
00318 if(clust[j] == trk.RecoPt(igap).GetPC()) infit[j] = 0;
00319 trk.Erase(igap);
00320 }
00321 } else {
00322 for (unsigned int i=0; i<igap; ++i) {
00323 for(unsigned int j=0; j<clust.size(); j++)
00324 if(clust[j] == trk.RecoPt(0).GetPC()) infit[j] = 0;
00325 trk.Erase(0);
00326 }
00327 }
00328 // Reset -> Check all gaps again after removel
00329 exgap=1.0E10;
00330 exhit=0;
00331 }
00332 }
00333
00334 if(trk.NRecoPt()<=1) break;
00335
00336 // 5. subtract hits from Hough map Subtract hits associated
00337 // with track for the next step to save calculation time.
00338 unsigned int icl = 0;
00339 for(unsigned int ih1=0; ih1<trk.NRecoPt()-1; ih1++) {
00340 for(unsigned int ih2=ih1+1; ih2<trk.NRecoPt(); ih2++) {
00341
00342 for(unsigned int j=0; j<clust.size(); j++)
00343 if(clust[j] == trk.RecoPt(ih1).GetPC()) icl = j;
00344
00345 hough.AddPoint(trk.RecoPt(ih1).GetPC()->Zpos()-Zave,
00346 trk.RecoPt(ih1).GetPC()->dZ(),
00347 trk.RecoPt(ih2).GetPC()->Zpos()-Zave,
00348 trk.RecoPt(ih2).GetPC()->dZ(),
00349 trk.RecoPt(ih1).GetPC()->Tpos()-Tave,
00350 trk.RecoPt(ih1).GetPC()->dT(),
00351 trk.RecoPt(ih2).GetPC()->Tpos()-Tave,
00352 trk.RecoPt(ih2).GetPC()->dT(),
00353 -1., icl);
00354 }
00355 }
00356
00357 if(trk.NRecoPt()<=1) break;
00358
00359 // Track candidates should have a minimum number of hits
00360 if(trk.NRecoPt()>=fMinTrackHits) {
00361 trk.SetXY(xychoice);
00362 track.push_back(trk);
00363 nfound++;
00364 }
00365 else {
00366 // Not enough hits, do not fill it.
00367 for(unsigned int j=0; j<clust.size(); j++)
00368 if(clust[j] == trk.RecoPt(0).GetPC()) infit[j] = 0;
00369 trk.Erase(0);
00370 }
00371
00372 if(!hough.UpdateTree()) break;
00373 }
00374
00375 return;
00376 }
|
|
|
Reimplemented from jobc::Module. Definition at line 57 of file FindTrackSeg.cxx. References geo::PlaneGeo::Cell(), DoView(), fDetGeom, fGeo, fInterval, geo::CellGeo::GetCenter(), edm::EventHandle::Header(), geo::Geometry::Instance(), geo::Geometry::Plane(), edm::EventHandle::Reco(), jobc::Result, and geo::PlaneGeo::View().
00058 {
00059
00060 // Get PlaneCluster
00061 vector<const recobase::PlaneCluster*> clust(0);
00062 try {
00063 evt.Reco().Get("./",clust);
00064 }
00065 catch (edm::Exception e) {
00066 return jobc::kFailed;
00067 }
00068
00069 // Prepare local variables
00070 if (!fGeo) fGeo = &geo::Geometry::Instance(evt.Header().Run(), fDetGeom);
00071
00072 double x1[3], x2[3];
00073 fGeo->Plane(0).Cell(0).GetCenter(x1);
00074 fGeo->Plane(2).Cell(0).GetCenter(x2);
00075 fInterval = x2[2] - x1[2];
00076
00077 vector<recobase::TrackBase> track(0);
00078
00079 vector<const recobase::PlaneCluster*> xclust(0);
00080 vector<const recobase::PlaneCluster*> yclust(0);
00081 for (unsigned int i=0; i<clust.size(); ++i) {
00082 if(fGeo->Plane(clust[i]->Plane()).View()==geo::kX) {
00083 xclust.push_back(clust[i]);
00084 } else {
00085 yclust.push_back(clust[i]);
00086 }
00087 }
00088
00089 this->DoView(xclust, track, geo::kX);
00090 this->DoView(yclust, track, geo::kY);
00091
00092 evt.Reco().MakeFolder("./Segment");
00093 for (unsigned int i=0; i<track.size(); ++i) {
00094 evt.Reco().Put(track[i],"./Segment");
00095 }
00096
00097 return jobc::kPassed;
00098 }
|
|
|
Implements cfg::Observer. Definition at line 46 of file FindTrackSeg.cxx. References fDetGeom, fGapThreshold, fHoughThreshold, fMaxTracks, and fMinTrackHits.
00047 {
00048 c("DetGeom"). Get(fDetGeom);
00049 c("HoughThreshold").Get(fHoughThreshold);
00050 c("MaxTracks"). Get(fMaxTracks);
00051 c("GapThreshold"). Get(fGapThreshold);
00052 c("MinTrackHits"). Get(fMinTrackHits);
00053 }
|
|
|
Definition at line 41 of file FindTrackSeg.h. |
|
|
Allowed gap sizes in tracks.
Definition at line 42 of file FindTrackSeg.h. |
|
|
Definition at line 46 of file FindTrackSeg.h. Referenced by Reco(). |
|
|
Threshold to apply to Hough map.
Definition at line 45 of file FindTrackSeg.h. |
|
|
Distance between planes.
Definition at line 47 of file FindTrackSeg.h. |
|
|
Max. tracks allowed per view.
Definition at line 43 of file FindTrackSeg.h. |
|
|
Min. hits allowed to form a track.
Definition at line 44 of file FindTrackSeg.h. |
1.3.5