#include <HoughView.h>
Public Member Functions | |
HoughView () | |
~HoughView () | |
virtual void | BuildDisplay (TCanvas *c) |
virtual void | DrawEvent (ParticleObjectHolder *poh=0, NtpStRecord *ntp=0, int ntpEvt=-1) |
void | DrawChains (int view) |
void | DrawTrueVertex (int view) |
void | DrawVertex (int view) |
void | DrawClusters (int view, TH2 *h) |
void | DrawHough (int view, TH2 *h) |
Private Member Functions | |
void | ClearHit (TH2D *his, double below_val, int curx, int cury) |
void | SaveHitGroup (TH2D *his, TH2D *saver, double save_val, double with_val, int curx, int cury) |
void | GetPeakAreaAverage (double &x, double &y, double &val, int &cnt, int curx, int cury, TH2D *hist) |
Private Attributes | |
TPad * | padU |
TPad * | padV |
TH2D * | viewU |
TH2D * | viewV |
TH2D * | houghviewU |
TH2D * | houghviewV |
ParticleObjectHolder * | mypoh |
Definition at line 13 of file HoughView.h.
HoughView::HoughView | ( | ) |
Definition at line 24 of file HoughView.cxx.
References houghviewU, houghviewV, padU, padV, viewU, and viewV.
00024 :Page() 00025 { 00026 viewU=0; 00027 viewV=0; 00028 padU=0; 00029 padV=0; 00030 houghviewU=0; 00031 houghviewV=0; 00032 00033 }
HoughView::~HoughView | ( | ) |
Definition at line 36 of file HoughView.cxx.
void HoughView::BuildDisplay | ( | TCanvas * | c | ) | [virtual] |
Reimplemented from Page.
Definition at line 40 of file HoughView.cxx.
References Page::myCanvas, padU, and padV.
Referenced by ParticleDisplay::BuildDisplay().
void HoughView::ClearHit | ( | TH2D * | his, | |
double | below_val, | |||
int | curx, | |||
int | cury | |||
) | [private] |
Definition at line 932 of file HoughView.cxx.
00933 { 00934 if(curx<1 || cury < 1 || curx> his->GetNbinsX() || cury >his->GetNbinsY())return; 00935 00936 // printf("looking in %d %d\n",curx,cury); 00937 00938 double thisval = his->GetBinContent(curx,cury); 00939 if(thisval>=below_val)return; 00940 00941 00942 // printf("val %f >= %f ?\n",thisval,below_val); 00943 ClearHit(his,below_val,curx,cury+1); 00944 ClearHit(his,below_val,curx+1,cury); 00945 ClearHit(his,below_val,curx,cury-1); 00946 ClearHit(his,below_val,curx-1,cury); 00947 ClearHit(his,below_val,curx+1,cury+1); 00948 ClearHit(his,below_val,curx-1,cury-1); 00949 ClearHit(his,below_val,curx+1,cury-1); 00950 ClearHit(his,below_val,curx-1,cury+1); 00951 00952 // printf("recursion done\n"); 00953 00954 his->SetBinContent(curx,cury,0); 00955 00956 }
void HoughView::DrawChains | ( | int | view | ) |
Definition at line 174 of file HoughView.cxx.
References ParticleObjectHolder::chu, ParticleObjectHolder::chv, copy(), Chain::end_t, Chain::end_z, ChainHelper::FindMaxPath(), ChainHelper::finished, ChainHelper::GetAllChildren(), ChainHelper::GetChain(), Chain::myId, mypoh, Chain::parentChain, Chain::t, and Chain::z.
00175 { 00176 00177 ChainHelper * ch = 0; 00178 if(view==2) 00179 ch=&mypoh->chu; 00180 else if(view==3) 00181 ch=&mypoh->chv; 00182 else return; 00183 00184 00185 std::vector<int> path =ch->FindMaxPath(); 00186 00187 for(unsigned int i=0;i<ch->finished.size();i++) 00188 { 00189 if(ch->finished[i].parentChain>-1 || ch->finished[i].entries < 1)continue; 00190 int parentcount=2; 00191 Chain *todraw=ch->GetChain(ch->finished[i].myId); 00192 00193 //draw the parent 00194 int dsize = todraw->t.size(); 00195 int needconnect=0; 00196 00197 if (dsize==1) //chain with only 1 cluster needs to draw line to parent explicitly... if not a head chain 00198 //other wise, make the single cluster, head chain, no children noticeable 00199 { 00200 // if(todraw->parentChain<0)continue; 00201 if(todraw->parentChain>-1) 00202 { 00203 dsize=2; 00204 needconnect=1; 00205 } 00206 } 00207 00208 double * tt = new double[todraw->t.size()]; 00209 double * tz = new double[todraw->t.size()]; 00210 if(needconnect) 00211 { 00212 tt[0]=ch->GetChain(todraw->parentChain)->end_t; 00213 tz[0]=ch->GetChain(todraw->parentChain)->end_z; 00214 tz[1]=todraw->z[0]; 00215 tt[1]=todraw->t[0]; 00216 } 00217 else 00218 { 00219 std::copy(todraw->t.begin(),todraw->t.end(),tt); 00220 std::copy(todraw->z.begin(),todraw->z.end(),tz); 00221 00222 if(dsize==1) 00223 { 00224 tt[1]=tt[0]; 00225 tz[1]=tz[0]-0.03; 00226 dsize=2; 00227 00228 } 00229 00230 } 00231 00232 TPolyLine *pl = new TPolyLine(dsize,tz,tt ); 00233 00234 //int linecol =2+todraw->level; 00235 00236 // pl->SetLineColor(linecol); 00237 // pl->SetLineColor(parentcount); 00238 pl->SetLineColor(2); 00239 00240 int width=1; 00241 //should make this better if it is to be really used 00242 00243 for(unsigned int j=0;j<path.size();j++) 00244 { 00245 if(path[j]==todraw->myId)width+=1; 00246 } 00247 00248 pl->SetLineWidth(width); 00249 // pl->Draw(); 00250 pl->Draw("same"); 00251 00252 for (unsigned int k=0;k<ch->finished[i].children.size();k++) 00253 { 00254 parentcount++; 00255 00256 std::vector<int> c = ch->GetAllChildren(ch->GetChain(ch->finished[i].children[k])); 00257 00258 printf("parent %d children ", ch->finished[i].myId); 00259 00260 for(unsigned int j=0;j<c.size();j++) 00261 { 00262 00263 Chain *todraw=ch->GetChain(c[j]); 00264 00265 printf("%d ",todraw->myId); 00266 00267 00268 int dsize = todraw->t.size(); 00269 int needconnect=0; 00270 00271 if (dsize==1) //chain with only 1 cluster needs to draw line to parent explicitly... if not a head chain 00272 //other wise, make the single cluster, head chain, no children noticeable 00273 { 00274 // if(todraw->parentChain<0)continue; 00275 00276 if(todraw->parentChain>-1) 00277 { 00278 dsize=2; 00279 needconnect=1; 00280 } 00281 } 00282 00283 double * tt = new double[todraw->t.size()]; 00284 double * tz = new double[todraw->t.size()]; 00285 if(needconnect) 00286 { 00287 tt[0]=ch->GetChain(todraw->parentChain)->end_t; 00288 tz[0]=ch->GetChain(todraw->parentChain)->end_z; 00289 tz[1]=todraw->z[0]; 00290 tt[1]=todraw->t[0]; 00291 } 00292 else 00293 { 00294 std::copy(todraw->t.begin(),todraw->t.end(),tt); 00295 std::copy(todraw->z.begin(),todraw->z.end(),tz); 00296 00297 if(dsize==1) 00298 { 00299 tt[1]=tt[0]; 00300 tz[1]=tz[0]-0.03; 00301 dsize=2; 00302 00303 } 00304 00305 } 00306 00307 TPolyLine *pl = new TPolyLine(dsize,tz,tt ); 00308 00309 //int linecol =2+todraw->level; 00310 00311 // pl->SetLineColor(linecol); 00312 pl->SetLineColor(parentcount); 00313 00314 00315 int width=1; 00316 //should make this better if it is to be really used 00317 00318 for(unsigned int j=0;j<path.size();j++) 00319 { 00320 if(path[j]==todraw->myId)width+=1; 00321 } 00322 00323 pl->SetLineWidth(width); 00324 // pl->Draw(); 00325 pl->Draw("same"); 00326 00327 }printf("\n"); 00328 }} 00329 00330 }
void HoughView::DrawClusters | ( | int | view, | |
TH2 * | h | |||
) |
Definition at line 336 of file HoughView.cxx.
References ParticleObjectHolder::cluster_saver, ParticleObjectHolder::clusterer, Managed::ManagedCluster::e, Managed::ClusterSaver::GetCluster(), Managed::ClusterManager::GetCluster(), Managed::ClusterSaver::GetClusterMap(), Managed::ClusterManager::GetClusterMap(), Managed::ManagedCluster::GetStatus(), Munits::m, mypoh, Managed::ManagedCluster::t, and Managed::ManagedCluster::z.
Referenced by DrawEvent().
00337 { 00338 00339 std::map<double, std::map<double, int> >::iterator p_iterr; 00340 std::map<double, int>::iterator s_iterr; 00341 00342 00343 Managed::ClusterManager *cluh = &mypoh->clusterer; 00344 /* 00345 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00346 { 00347 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00348 { 00349 00350 00351 h->Fill( p_iterr->first,s_iterr->first, cluh->GetCluster(s_iterr->second)->e); 00352 } 00353 } 00354 h->SetStats(0); 00355 h->Draw("colz"); 00356 */ 00357 00358 //Managed::ManagedCluster *largest=0; 00359 //double last_plane=0; 00360 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00361 { 00362 00363 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00364 { 00365 //Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00366 h->Fill( p_iterr->first,s_iterr->first, cluh->GetCluster(s_iterr->second)->e); 00367 00368 } 00369 } 00370 00371 00372 00373 Managed::ClusterSaver * cluhs = &mypoh->cluster_saver; 00374 for(p_iterr=cluhs->GetClusterMap(view)->begin();p_iterr!=cluhs->GetClusterMap(view)->end(); p_iterr++) 00375 { 00376 00377 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00378 { 00379 Managed::ManagedCluster * clus = cluhs->GetCluster(s_iterr->second); 00380 00381 //printf("%f %f %d\n",clus->z,clus->t,clus->GetStatus()); 00382 if(clus->GetStatus()==1)continue;//long muon clusters 00383 00384 h->Fill( p_iterr->first,s_iterr->first, cluhs->GetCluster(s_iterr->second)->e); 00385 00386 } 00387 } 00388 h->SetStats(0); 00389 h->Draw("colz"); 00390 00391 00392 00393 /* 00394 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00395 { 00396 //new plane and a largest cluster in the previous plane? 00397 if(fabs(last_plane-p_iterr->first)>0.03 && largest) 00398 { 00399 TMarker *m=new TMarker(largest->z,largest->t, 25); 00400 m->Draw("same"); 00401 largest=0; 00402 last_plane=p_iterr->first; 00403 } 00404 00405 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00406 { 00407 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00408 00409 if(!largest) 00410 { 00411 largest=clus; 00412 continue; 00413 } 00414 if(largest->e<clus->e)largest=clus; 00415 } 00416 } 00417 if(fabs(last_plane-largest->z)>0.03 && largest) 00418 { 00419 TMarker *m=new TMarker(largest->z,largest->t, 25); 00420 m->Draw("same"); 00421 largest=0; 00422 } 00423 */ 00424 00425 //to take all within xxx% of the max hit in the event! 00426 00427 //find max value! 00428 double maxvale=0; 00429 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00430 { 00431 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00432 { 00433 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00434 if(clus->e>maxvale) 00435 { 00436 maxvale=clus->e; 00437 } 00438 } 00439 } 00440 00441 00442 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00443 { 00444 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00445 { 00446 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00447 00448 // if(clus->e>maxvale*0.1) 00449 // { 00450 if(clus->GetStatus()==1)continue;//long muon clusters 00451 TMarker *m=new TMarker(clus->z,clus->t, 25); 00452 m->Draw("same"); 00453 00454 // } 00455 } 00456 } 00457 00458 00459 /* 00460 for(p_iterr=cluhs->GetClusterMap(view)->begin();p_iterr!=cluhs->GetClusterMap(view)->end(); p_iterr++) 00461 { 00462 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00463 { 00464 Managed::ManagedCluster * clus = cluhs->GetCluster(s_iterr->second); 00465 00466 if(clus->GetStatus()==1)continue;//long muon clusters 00467 TMarker *m=new TMarker(clus->z,clus->t, 27); 00468 m->Draw("same"); 00469 00470 00471 } 00472 } 00473 */ 00474 00475 00476 00477 }
void HoughView::DrawEvent | ( | ParticleObjectHolder * | poh = 0 , |
|
NtpStRecord * | ntp = 0 , |
|||
int | ntpEvt = -1 | |||
) | [virtual] |
Reimplemented from Page.
Definition at line 56 of file HoughView.cxx.
References DrawClusters(), DrawHough(), DrawTrueVertex(), DrawVertex(), ParticleObjectHolder::event, RecRecordImp< T >::GetHeader(), VldContext::GetSimFlag(), RecHeader::GetVldContext(), houghviewU, houghviewV, isMC, SimFlag::kMC, ParticleEvent::maxu, ParticleEvent::maxv, ParticleEvent::maxz, ParticleObjectHolder::mctrue, ParticleEvent::minu, ParticleEvent::minv, ParticleEvent::minz, Page::myCanvas, mypoh, padU, padV, viewU, viewV, MCTrue::vtx_u, MCTrue::vtx_v, and MCTrue::vtx_z.
Referenced by ParticleDisplay::UpdateDisplay().
00057 { 00058 myCanvas->cd(); 00059 00060 padU->Clear(); 00061 padV->Clear(); 00062 00063 padU->Divide(2,1); 00064 padV->Divide(2,1); 00065 00066 mypoh=poh; 00067 00068 int isMC=0; 00069 if(poh->GetHeader().GetVldContext().GetSimFlag() == SimFlag::kMC)isMC=1; 00070 00071 double minz= mypoh->event.minz; 00072 double minu= mypoh->event.minu; 00073 double minv= mypoh->event.minv; 00074 double maxz= mypoh->event.maxz; 00075 double maxu= mypoh->event.maxu; 00076 double maxv= mypoh->event.maxv; 00077 00078 00079 00080 if(isMC) 00081 { 00082 if(mypoh->mctrue.vtx_z>=0)minz=minz<mypoh->mctrue.vtx_z ? minz : mypoh->mctrue.vtx_z; 00083 maxz=maxz>mypoh->mctrue.vtx_z ? maxz : mypoh->mctrue.vtx_z; 00084 minu=minu<mypoh->mctrue.vtx_u ? minu : mypoh->mctrue.vtx_u; 00085 maxu=maxu>mypoh->mctrue.vtx_u ? maxu : mypoh->mctrue.vtx_u; 00086 minv=minv<mypoh->mctrue.vtx_v ? minv : mypoh->mctrue.vtx_v; 00087 maxv=maxv>mypoh->mctrue.vtx_v ? maxv : mypoh->mctrue.vtx_v; 00088 } 00089 00090 00091 00092 00093 minz -= 0.1; 00094 minu -= 0.1; 00095 minv -= 0.1; 00096 maxz += 0.1; 00097 maxu += 0.1; 00098 maxv += 0.1; 00099 00100 int nu=(int)((maxu-minu)/0.0412); 00101 int nv=(int)((maxv-minv)/0.0412); 00102 int nz=(int)((maxz-minz)/0.0354); 00103 00104 00105 if(viewU)delete viewU; 00106 if(viewV)delete viewV; 00107 00108 if(houghviewU)delete houghviewU; 00109 if(houghviewV)delete houghviewV; 00110 00111 viewU = new TH2D("hu","Clusters U", nz, minz,maxz,nu,minu, maxu); 00112 viewV = new TH2D("hv","Clusters V", nz, minz,maxz,nv,minv, maxv); 00113 houghviewU = new TH2D("hhu","Hough U", 80, 0,3.141592,50,-1,1); 00114 houghviewV = new TH2D("hhv","Hough V", 80, 0,3.141592,50,-1,1); 00115 viewU->SetStats(0); 00116 viewV->SetStats(0); 00117 viewU->SetDirectory(0); 00118 viewV->SetDirectory(0); 00119 00120 houghviewU->SetStats(0); 00121 houghviewV->SetStats(0); 00122 houghviewU->SetDirectory(0); 00123 houghviewV->SetDirectory(0); 00124 00125 if(mypoh) 00126 { 00127 00128 00129 //psf.FindPrimaryShower(&mypoh->clusterer); 00130 00131 00132 padU->cd(1); 00133 00134 viewU->Draw(); 00135 DrawClusters(2,viewU); 00136 DrawVertex(2); 00137 if(isMC)DrawTrueVertex(2); 00138 // DrawChains(2); 00139 00140 padU->cd(2); 00141 DrawHough(2,houghviewU); 00142 00143 00144 00145 padV->cd(1); 00146 viewV->Draw(); 00147 DrawClusters(3,viewV); 00148 DrawVertex(3); 00149 if(isMC)DrawTrueVertex(3); 00150 // DrawChains(3); 00151 padV->cd(2); 00152 DrawHough(3,houghviewV); 00153 00154 00155 00156 00157 } 00158 00159 00160 double maxz_u = viewU->GetMaximum(); 00161 double maxz_v = viewV->GetMaximum(); 00162 double maxz_both = maxz_u>maxz_v?maxz_u:maxz_v; 00163 viewU->GetZaxis()->SetRangeUser(0, maxz_both); 00164 viewV->GetZaxis()->SetRangeUser(0, maxz_both); 00165 00166 padU->Update(); 00167 padV->Update(); 00168 00169 myCanvas->Update(); 00170 00171 }
void HoughView::DrawHough | ( | int | view, | |
TH2 * | h | |||
) |
the old code....
Definition at line 518 of file HoughView.cxx.
References PrimaryShowerFinder::FindFirstIntersection(), PrimaryShowerFinder::GetHoughLines(), PrimaryShowerFinder::intU, PrimaryShowerFinder::intV, Munits::m, max, mypoh, padU, padV, ParticleObjectHolder::psf, and PrimaryShowerFinder::ran.
Referenced by DrawEvent().
00519 { 00520 00521 if(view==2)padU->cd(2); 00522 if(view==3)padV->cd(2); 00523 00524 double t=0; 00525 double z=0; 00526 00527 if(!mypoh->psf.ran)return; 00528 int cp = mypoh->psf.FindFirstIntersection(view,t,z); 00529 if(cp) 00530 { 00531 TMarker *m=new TMarker(z,t, 29); 00532 m->Draw("same"); 00533 printf("closest intersection z %f t %f\n",z,t); 00534 } 00535 00536 00537 00538 TH2D * inth = view==2? mypoh->psf.intU:mypoh->psf.intV; 00539 if(inth) 00540 { 00541 inth->SetStats(0); 00542 inth->Draw("colz"); 00543 } 00544 00545 /* 00546 00547 TH2D * hh = mypoh->psf.GetHoughMap(view); 00548 hh->SetStats(0); 00549 hh->Draw("colz"); 00550 */ 00551 00552 if(view==2)padU->cd(1); 00553 if(view==3)padV->cd(1); 00554 00555 00556 00557 00558 //double view_xmin = view==2?viewU->GetXaxis()->GetXmin():viewV->GetXaxis()->GetXmin(); 00559 //double view_xmax = view==2?viewU->GetXaxis()->GetXmax():viewV->GetXaxis()->GetXmax(); 00560 00561 00562 00563 00564 00565 00566 std::vector<HoughLine>*houghlines = mypoh->psf.GetHoughLines(view); 00567 00568 printf("# hough lines %d\n",(int)houghlines->size()); 00569 00570 int max =-1; 00571 if((*houghlines).size()>50)max = (*houghlines).size()-50; 00572 for(int i=(*houghlines).size()-1;i>max;i--) 00573 { 00574 if((*houghlines)[i].ncluster<1)continue; 00575 // printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f offz %f expt %f phi %f\n", (*houghlines)[i].sum_e, (*houghlines)[i].ncluster, (*houghlines)[i].start_t, (*houghlines)[i].start_z, (*houghlines)[i].end_t, (*houghlines)[i].end_z, (*houghlines)[i].chi2/(*houghlines)[i].ncluster,(*houghlines)[i].offset_z,(*houghlines)[i].GetExpectedT((*houghlines)[i].offset_z),(*houghlines)[i].phi); 00576 00577 /* 00578 double ym = (*houghlines)[i].GetExpectedT(view_xmin); 00579 double yp = (*houghlines)[i].GetExpectedT(view_xmax); 00580 00581 00582 TLine *l=new TLine(view_xmin,ym,view_xmax,yp); 00583 */ 00584 00585 TLine *l=new TLine((*houghlines)[i].start_z,(*houghlines)[i].start_t,(*houghlines)[i].end_z,(*houghlines)[i].end_t); 00586 00587 if((*houghlines)[i].primary)l->SetLineWidth(2); 00588 l->Draw(); 00589 00590 // printf("fit line --%f %f - %f %f -- theta r %f %f\n",view_xmin,ym,view_xmax,yp,houghlines[i].theta,houghlines[i].r); 00591 00592 } 00593 00594 00595 00596 00597 00598 00599 00600 00601 00602 00603 00604 return; 00605 00607 /* 00608 00609 std::map<double, std::map<double, int> >::iterator p_iterr; 00610 std::map<double, int>::iterator s_iterr; 00611 00612 printf("----view %d----\n\n",view); 00613 Managed::ClusterManager *cluh = view==2?&mypoh->clusterer:&mypoh->clusterer; 00614 00615 00616 double off_z = cluh->GetClusterMap(view)->begin()->first; 00617 double off_t = cluh->GetClusterMap(view)->begin()->second.begin()->first; 00618 00619 00620 / * 00621 for(p_iterr=cluh->cluster_map.begin();p_iterr!=cluh->cluster_map.end(); p_iterr++) 00622 { 00623 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00624 { 00625 00626 for(double i=0;i<2*3.141592;i+=3.141592*2/99) 00627 { 00628 h->Fill(i, (p_iterr->first-off_z) * cos(i) + (s_iterr->first-off_t) * sin(i), 1cluh->GetCluster(s_iterr->second)->e); 00629 } 00630 } 00631 } 00632 00633 *//* 00634 00635 Managed::ManagedCluster *largest=0; 00636 double last_plane=0; 00637 00638 //to take the largest in each plane 00639 / * 00640 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00641 { 00642 //new plane and a largest cluster in the previous plane? 00643 if(fabs(last_plane-p_iterr->first)>0.03 && largest) 00644 { 00645 for(double i=0;i<3.141592;i+=3.141592/(h->GetNbinsX()-1)) 00646 { 00647 //double val = largest->e; 00648 double val=1; 00649 h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i), val); 00650 00651 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+2./(h->GetNbinsY())., val*.5); 00652 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+4./(h->GetNbinsY())., val*.5*.5); 00653 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-2./(h->GetNbinsY())., val*.5); 00654 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-4./(h->GetNbinsY())., val*.5*.5); 00655 } 00656 largest=0; 00657 last_plane=p_iterr->first; 00658 } 00659 00660 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00661 { 00662 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00663 if(!largest) 00664 { 00665 largest=clus; 00666 continue; 00667 } 00668 if(largest->e<clus->e)largest=clus; 00669 } 00670 } 00671 if(fabs(last_plane-largest->z)>0.03 && largest) 00672 { 00673 for(double i=0;i<3.141592;i+=3.141592/(h->GetNbinsX()-1)) 00674 { 00675 //double val = largest->e; 00676 double val=1; 00677 h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i), val); 00678 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+2./(h->GetNbinsY())., val*.5); 00679 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+4./(h->GetNbinsY())., val*.5*.5); 00680 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-2./(h->GetNbinsY())., val*.5); 00681 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-4./(h->GetNbinsY())., val*.5*.5); 00682 00683 } 00684 00685 largest=0; 00686 } 00687 */ /* 00688 00689 00690 TH2D copy; 00691 h->Copy(copy); 00692 00693 //to take all within xxx% of the max hit in the event! 00694 00695 //find max value! 00696 double maxvale=0; 00697 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00698 { 00699 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00700 { 00701 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00702 if(clus->e>maxvale) 00703 { 00704 maxvale=clus->e; 00705 } 00706 } 00707 } 00708 00709 00710 for(p_iterr=cluh->GetClusterMap(view)->begin();p_iterr!=cluh->GetClusterMap(view)->end(); p_iterr++) 00711 { 00712 for(s_iterr=p_iterr->second.begin();s_iterr!=p_iterr->second.end(); s_iterr++) 00713 { 00714 Managed::ManagedCluster * clus = cluh->GetCluster(s_iterr->second); 00715 00716 // if(clus->e>maxvale*0.1) 00717 // { 00718 for(double i=0;i<3.141592;i+=3.141592/(h->GetNbinsX()-1)) 00719 { 00720 //double val = clus->e; 00721 double val=1; 00722 h->Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i), val); 00723 copy.Fill(i, (clus->z-off_z) * cos(i) + (clus->t-off_t) * sin(i), clus->e); 00724 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+2./(h->GetNbinsY())., val*.5); 00725 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)+4./(h->GetNbinsY())., val*.5*.5); 00726 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-2./(h->GetNbinsY())., val*.5); 00727 // h->Fill(i, (largest->z-off_z) * cos(i) + (largest->t-off_t) * sin(i)-4./(h->GetNbinsY())., val*.5*.5); 00728 00729 } 00730 // } 00731 } 00732 } 00733 00734 00735 / * 00736 h->SetStats(0); 00737 h->Draw("colz"); 00738 00739 return; 00740 */ 00741 //peak filter the hough map! 00742 /* 00743 TH2D ch; 00744 h->Copy(ch); 00745 00746 00747 h->Reset(); 00748 00749 int nx=ch.GetNbinsX(); 00750 int ny=ch.GetNbinsY(); 00751 00752 while(ch.GetMaximum()>1) 00753 { 00754 int x; 00755 int y; 00756 int z; 00757 int m = ch.GetMaximumBin(); 00758 ch.GetBinXYZ(m,x,y,z); 00759 00760 / * printf("clearing hits around %d %d\n",x,y); 00761 ClearHit(&ch, ch.GetMaximum(), x+1,y); 00762 ClearHit(&ch, ch.GetMaximum(), x,y+1); 00763 ClearHit(&ch, ch.GetMaximum(), x,y-1); 00764 ClearHit(&ch, ch.GetMaximum(), x-1,y); 00765 ClearHit(&ch, ch.GetMaximum(), x+1,y+1); 00766 ClearHit(&ch, ch.GetMaximum(), x-1,y-1); 00767 ClearHit(&ch, ch.GetMaximum(), x+1,y-1); 00768 ClearHit(&ch, ch.GetMaximum(), x-1,y+1); 00769 */ /* 00770 SaveHitGroup(&ch,(TH2D*)h, (double)ch.GetBinContent(x,y),(double)ch.GetBinContent(x,y), x,y); 00771 00772 00773 } 00774 00775 h->SetStats(0); 00776 h->Draw("colz"); 00777 00778 00779 if(view==2)padU->cd(1); 00780 if(view==3)padV->cd(1); 00781 00782 00783 / * 00784 for(int i=1;i<=h->GetNbinsX();i++) 00785 for(int j=1;j<=h->GetNbinsY();j++) 00786 { 00787 if(h->GetBinContent(i,j)>=4) 00788 { 00789 00790 double theta = h->GetXaxis()->GetBinLowEdge(i); 00791 double r = h->GetYaxis()->GetBinLowEdge(j); 00792 00793 double ym = (- cos(theta)/sin(theta))*(viewU->GetXaxis()->GetXmin()-off_z)+r/sin(theta) + off_t; 00794 double yp = (- cos(theta)/sin(theta))*(viewU->GetXaxis()->GetXmax()-off_z)+r/sin(theta) + off_t; 00795 00796 TLine *l=new TLine(viewU->GetXaxis()->GetXmin(),ym,viewU->GetXaxis()->GetXmax(),yp); 00797 l->Draw(); 00798 00799 } 00800 } 00801 */ 00802 /* 00803 multimap<double, std::pair<double,double> > top_map; 00804 00805 int tops=0; 00806 / * for(int i=1;i<=h->GetNbinsX();i++) 00807 for(int j=1;j<=h->GetNbinsY();j++) 00808 { 00809 if(h->GetBinContent(i,j)>1) 00810 { 00811 top_map.insert(make_pair((double)h->GetBinContent(i,j), std::pair<int,int>(i,j) )); 00812 tops++; 00813 printf("topbin found x %d y %d\n",i,j); 00814 } 00815 } 00816 */ 00817 00818 //h->Rebin2D(4,4); 00819 /* while(h->GetMaximum()>0) 00820 { 00821 double xv=0; 00822 double yv=0; 00823 double val=0; 00824 int cnt =0; 00825 00826 int x; 00827 int y; 00828 int z; 00829 int m = h->GetMaximumBin(); 00830 h->GetBinXYZ(m,x,y,z); 00831 00832 GetPeakAreaAverage(xv, yv,val, cnt, x, y, (TH2D*)h); 00833 xv/=(double)cnt; 00834 yv/=(double)cnt; 00835 00836 val=copy.GetBinContent(copy.FindBin(xv,yv)); 00837 00838 top_map.insert(make_pair((double)val, std::pair<double,double>(xv,yv)) ); 00839 tops++; 00840 // printf("topbin found x %f y %f\n",xv,yv); 00841 00842 00843 } 00844 00845 00846 00847 printf("printing top %d\n",tops); 00848 00849 multimap<double, std::pair<double,double> >::reverse_iterator it = top_map.rbegin(); 00850 00851 double view_xmin = view==2?viewU->GetXaxis()->GetXmin():viewV->GetXaxis()->GetXmin(); 00852 double view_xmax = view==2?viewU->GetXaxis()->GetXmax():viewV->GetXaxis()->GetXmax(); 00853 00854 std::vector<HoughLine> houghlines; 00855 00856 for(unsigned int k=0;k<tops;k++) 00857 { 00858 //int i=it->second.first; 00859 //int j=it->second.second; 00860 00861 //double theta = h->GetXaxis()->GetBinCenter(i); 00862 //double r = h->GetYaxis()->GetBinCenter(j); 00863 double theta=it->second.first; 00864 double r = it->second.second; 00865 00866 // double ym = (- cos(theta)/sin(theta))*(view_xmin-off_z)+r/sin(theta) + off_t; 00867 // double yp = (- cos(theta)/sin(theta))*(view_xmax-off_z)+r/sin(theta) + off_t; 00868 00869 00870 HoughLine hl(theta, r, off_t, off_z); 00871 00872 PrimaryShowerFinder::LoadCloseHits(&hl,&mypoh->clusterer,view,0.03); 00873 00874 houghlines.push_back(hl); 00875 00876 it++; 00877 } 00878 00879 sort(&houghlines[0],&houghlines[houghlines.size()],CompareLength); 00880 00881 //should remove duplicates --- (caused by the same clusters being matched to different hough map peaks) 00882 std::vector<HoughLine> houghlinestemp; 00883 for(int i=0;i<houghlines.size();i++) 00884 { 00885 HoughLine * l = &houghlines[i]; 00886 HoughLine * r=0; 00887 if(i<houghlines.size()-1)r = &houghlines[i+1]; 00888 if(houghlines[i].ncluster>1)houghlinestemp.push_back(houghlines[i]); 00889 if(l && r && l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t) 00890 { 00891 00892 while (i<houghlines.size()) 00893 { 00894 i++; 00895 HoughLine * l = &houghlines[i]; 00896 HoughLine * r = &houghlines[i+1]; 00897 if(! (l->start_z == r->start_z && l->start_t == r->start_t && l ->end_z == r->end_z && l->end_t == r->end_t)) break; 00898 } 00899 } 00900 00901 } 00902 00903 houghlines = houghlinestemp; 00904 00905 sort(&houghlines[0],&houghlines[houghlines.size()],CompareForwardAndClusters); 00906 00907 int max =-1; 00908 if(houghlines.size()>10)max = houghlines.size()-10; 00909 for(int i=houghlines.size()-1;i>max;i--) 00910 { 00911 if(houghlines[i].ncluster<1)break; 00912 printf("houghline with e %f clusters %d (t,z) (%f,%f) (%f,%f) chi2/ndf %f\n", houghlines[i].sum_e, houghlines[i].ncluster, houghlines[i].start_t, houghlines[i].start_z, houghlines[i].end_t, houghlines[i].end_z, houghlines[i].chi2/houghlines[i].ncluster); 00913 00914 double ym = houghlines[i].GetExpectedT(view_xmin); 00915 double yp = houghlines[i].GetExpectedT(view_xmax); 00916 00917 00918 TLine *l=new TLine(view_xmin,ym,view_xmax,yp); 00919 00920 if(i==houghlines.size()-1)l->SetLineWidth(2); 00921 l->Draw(); 00922 00923 // printf("fit line --%f %f - %f %f -- theta r %f %f\n",view_xmin,ym,view_xmax,yp,houghlines[i].theta,houghlines[i].r); 00924 00925 } 00926 00927 00928 */ 00929 }
void HoughView::DrawTrueVertex | ( | int | view | ) |
Definition at line 502 of file HoughView.cxx.
References Munits::m, ParticleObjectHolder::mctrue, mypoh, MCTrue::vtx_u, MCTrue::vtx_v, and MCTrue::vtx_z.
Referenced by DrawEvent().
void HoughView::DrawVertex | ( | int | view | ) |
Definition at line 487 of file HoughView.cxx.
References ParticleObjectHolder::event, Munits::m, mypoh, ParticleEvent::vtx_u, ParticleEvent::vtx_v, and ParticleEvent::vtx_z.
Referenced by DrawEvent().
void HoughView::GetPeakAreaAverage | ( | double & | x, | |
double & | y, | |||
double & | val, | |||
int & | cnt, | |||
int | curx, | |||
int | cury, | |||
TH2D * | hist | |||
) | [private] |
Definition at line 994 of file HoughView.cxx.
00995 { 00996 00997 double thisval = hist->GetBinContent(curx,cury); 00998 00999 if(thisval==0)return; 01000 01001 val+=thisval; 01002 x+=hist->GetXaxis()->GetBinCenter(curx); 01003 y+=hist->GetYaxis()->GetBinCenter(cury); 01004 cnt++; 01005 hist->SetBinContent(curx,cury,0); 01006 01007 01008 GetPeakAreaAverage(x, y, val,cnt, curx+1, cury, hist); 01009 GetPeakAreaAverage(x, y, val,cnt, curx, cury+1, hist); 01010 GetPeakAreaAverage(x, y, val,cnt, curx, cury-1, hist); 01011 GetPeakAreaAverage(x, y, val,cnt, curx-1, cury, hist); 01012 GetPeakAreaAverage(x, y, val,cnt, curx+1, cury+1, hist); 01013 GetPeakAreaAverage(x, y, val,cnt, curx-1, cury-1, hist); 01014 GetPeakAreaAverage(x, y, val,cnt, curx+1, cury-1, hist); 01015 GetPeakAreaAverage(x, y, val,cnt, curx-1, cury+1, hist); 01016 01017 01018 }
void HoughView::SaveHitGroup | ( | TH2D * | his, | |
TH2D * | saver, | |||
double | save_val, | |||
double | with_val, | |||
int | curx, | |||
int | cury | |||
) | [private] |
Definition at line 958 of file HoughView.cxx.
00959 { 00960 if(curx<1 || cury < 1 || curx> his->GetNbinsX() || cury >his->GetNbinsY())return; 00961 00962 00963 00964 double thisval = his->GetBinContent(curx,cury); 00965 00966 // printf("saving %d %d thisval %f saveval %f withval %f\n",curx,cury,thisval,save_val, with_val); 00967 if(thisval==0)return; 00968 if(fabs(thisval-save_val)<0.001) 00969 { 00970 saver->SetBinContent(curx, cury,thisval); 00971 //thisval=0; 00972 } 00973 else if(thisval>=with_val)return; 00974 00975 00976 his->SetBinContent(curx,cury,0); 00977 00978 SaveHitGroup(his,saver,save_val,thisval,curx+1,cury); 00979 SaveHitGroup(his,saver,save_val,thisval,curx,cury+1); 00980 SaveHitGroup(his,saver,save_val,thisval,curx-1,cury); 00981 SaveHitGroup(his,saver,save_val,thisval,curx,cury-1); 00982 SaveHitGroup(his,saver,save_val,thisval,curx+1,cury+1); 00983 SaveHitGroup(his,saver,save_val,thisval,curx-1,cury-1); 00984 SaveHitGroup(his,saver,save_val,thisval,curx+1,cury-1); 00985 SaveHitGroup(his,saver,save_val,thisval,curx-1,cury+1); 00986 // printf("save recursion done\n"); 00987 00988 00989 //printf("saving points %d %d %f\n",curx,cury,thisval); 00990 00991 }
TH2D* HoughView::houghviewU [private] |
Definition at line 39 of file HoughView.h.
Referenced by DrawEvent(), and HoughView().
TH2D* HoughView::houghviewV [private] |
Definition at line 40 of file HoughView.h.
Referenced by DrawEvent(), and HoughView().
ParticleObjectHolder* HoughView::mypoh [private] |
Definition at line 42 of file HoughView.h.
Referenced by DrawChains(), DrawClusters(), DrawEvent(), DrawHough(), DrawTrueVertex(), and DrawVertex().
TPad* HoughView::padU [private] |
Definition at line 32 of file HoughView.h.
Referenced by BuildDisplay(), DrawEvent(), DrawHough(), and HoughView().
TPad* HoughView::padV [private] |
Definition at line 33 of file HoughView.h.
Referenced by BuildDisplay(), DrawEvent(), DrawHough(), and HoughView().
TH2D* HoughView::viewU [private] |
Definition at line 36 of file HoughView.h.
Referenced by DrawEvent(), and HoughView().
TH2D* HoughView::viewV [private] |
Definition at line 37 of file HoughView.h.
Referenced by DrawEvent(), and HoughView().