WCSimRingFinder.cc

Go to the documentation of this file.
00001 #include "WCSimRingFinder.hh"
00002 
00003 #include "WCSimRecoDigit.hh"
00004 #include "WCSimRecoEvent.hh"
00005 #include "WCSimRecoVertex.hh"
00006 #include "WCSimRecoRing.hh"
00007 #include "WCSimHoughTransform.hh"
00008 #include "WCSimHoughTransformArray.hh"
00009 
00010 #include "WCSimGeometry.hh"
00011 
00012 #include <cmath>
00013 #include <iostream>
00014 #include <cassert>
00015 
00016 ClassImp(WCSimRingFinder)
00017 
00018 static WCSimRingFinder* fgRingFinder = 0;
00019 
00020 WCSimRingFinder* WCSimRingFinder::Instance()
00021 {
00022   if( !fgRingFinder ){
00023     fgRingFinder = new WCSimRingFinder();
00024   }
00025 
00026   if( !fgRingFinder ){
00027     assert(fgRingFinder);
00028   }
00029 
00030   if( fgRingFinder ){
00031 
00032   }
00033 
00034   return fgRingFinder;
00035 }
00036   
00037 void WCSimRingFinder::UseRecoVertex()
00038 {
00039   WCSimRingFinder::Instance()->SetUsingRecoVertex();
00040 }
00041 
00042 void WCSimRingFinder::HoughX( Int_t x )
00043 {
00044   WCSimRingFinder::Instance()->SetHoughX(x);
00045 }
00046   
00047 void WCSimRingFinder::HoughY( Int_t y )
00048 {
00049   WCSimRingFinder::Instance()->SetHoughY(y);
00050 }
00051   
00052 void WCSimRingFinder::HoughPoints( Int_t n )
00053 {
00054   WCSimRingFinder::Instance()->SetHoughPoints(n);
00055 }
00056   
00057 void WCSimRingFinder::ConeAngleBins( Int_t bins )
00058 {
00059   WCSimRingFinder::Instance()->SetConeAngleBins(bins);
00060 }
00061 
00062 void WCSimRingFinder::ConeAngleMin( Double_t min )
00063 {
00064   WCSimRingFinder::Instance()->SetConeAngleMin(min);
00065 }
00066 
00067 void WCSimRingFinder::ConeAngleMax( Double_t max )
00068 {
00069   WCSimRingFinder::Instance()->SetConeAngleMax(max);
00070 }
00071 
00072 void WCSimRingFinder::PrintParameters()
00073 {
00074   WCSimRingFinder::Instance()->RunPrintParameters();
00075 }
00076 
00077 WCSimRingFinder::WCSimRingFinder()
00078 {
00079   // use reco vertex
00080   fUseVertex = 0;
00081 
00082   // default hough transform parameters
00083   fHoughX = 120;          // bins in phi
00084   fHoughY = 60;           // bins in cos(theta)
00085   fHoughPoints = 360;     // hough points
00086 
00087   fConeAngleBins = 30;   // hough angle bins
00088   fConeAngleMin = 24.5;  // hough cone (min)
00089   fConeAngleMax = 54.5;  // hough cone (max)
00090   
00091   // hough transform object
00092   fHoughTransform = 0;
00093 
00094   // hough transform array
00095   fHoughTransformArray = 0;
00096 
00097   // make list of rings
00098   fRingList = new std::vector<WCSimRecoRing*>;
00099 }
00100 
00101 WCSimRingFinder::~WCSimRingFinder()
00102 {
00103   if( fHoughTransform ){
00104     delete fHoughTransform;
00105   }
00106 
00107   if( fHoughTransformArray ){
00108     delete fHoughTransformArray;
00109   }
00110 
00111   for( UInt_t i=0; i<fRingList->size(); i++ ){
00112     delete (WCSimRecoRing*)(fRingList->at(i));
00113   }
00114  
00115   fRingList->clear();
00116 
00117   delete fRingList;
00118 }
00119 
00120 void WCSimRingFinder::RunPrintParameters()
00121 {
00122   std::cout << " *** WCSimRingFinder::PrintParameters() *** " << std::endl;
00123 
00124   std::cout << "  Ring Finder Parameters: " << std::endl
00125             << "   UseVertex = " << fUseVertex << std::endl
00126             << "   HoughX = " << fHoughX << std::endl
00127             << "   HoughY = " << fHoughY << std::endl
00128             << "   ConeAngleBins = " << fConeAngleBins << std::endl
00129             << "   ConeAngleMin = " << fConeAngleMin << std::endl
00130             << "   ConeAngleMax = " << fConeAngleMax << std::endl;
00131 
00132   return;
00133 }
00134 
00135 void WCSimRingFinder::Reset()
00136 {
00137 
00138   return;
00139 }
00140 
00141 std::vector<WCSimRecoRing*>* WCSimRingFinder::Run(WCSimRecoVertex* myVertex)
00142 {
00143   std::cout << " *** WCSimRingFinder::Run(Vertex) *** " << std::endl;
00144 
00145   // reset ring list
00146   // ===============  
00147   for( UInt_t i=0; i<fRingList->size(); i++ ){
00148     delete (WCSimRecoRing*)(fRingList->at(i));
00149   }
00150  
00151   // clear ring list
00152   // ===============
00153   fRingList->clear();  
00154 
00155   // Make new ring, using vertex
00156   // ===========================
00157   WCSimRecoRing* myRing = new WCSimRecoRing(myVertex->GetX(),
00158                                             myVertex->GetY(),
00159                                             myVertex->GetZ(),   
00160                                             myVertex->GetDirX(),
00161                                             myVertex->GetDirY(),
00162                                             myVertex->GetDirZ(),
00163                                             myVertex->GetConeAngle(), 
00164                                             0.0); // height of hough peak
00165 
00166   fRingList->push_back(myRing);
00167 
00168   // Return Ring List
00169   // ================
00170   return fRingList;
00171 }
00172 
00173 std::vector<WCSimRecoRing*>* WCSimRingFinder::Run(WCSimRecoEvent* myEvent, WCSimRecoVertex* myVertex)
00174 {
00175   // get filter digit list
00176   // =====================
00177   std::vector<WCSimRecoDigit*>* myFilterDigitList = (std::vector<WCSimRecoDigit*>*)(myEvent->GetFilterDigitList());
00178 
00179   // run vertex finder
00180   // =================
00181   return (std::vector<WCSimRecoRing*>*)(this->Run(myFilterDigitList,myVertex));
00182 }
00183 
00184 std::vector<WCSimRecoRing*>* WCSimRingFinder::Run(std::vector<WCSimRecoDigit*>* myDigitList, WCSimRecoVertex* myVertex)
00185 {
00186   std::cout << " *** WCSimRingFinder::Run(...) *** " << std::endl;
00187 
00188   // override 
00189   // ========
00190   if( fUseVertex ){
00191     std::cout << "  --- reconstruct ring from vertex --- " << std::endl;
00192     return this->Run(myVertex);
00193   }
00194 
00195   // reset ring list
00196   // ===============  
00197   for( UInt_t i=0; i<fRingList->size(); i++ ){
00198     delete (WCSimRecoRing*)(fRingList->at(i));
00199   }
00200  
00201   // clear ring list
00202   // ===============
00203   fRingList->clear();
00204 
00205   // apply Hough Transform
00206   // =====================
00207   WCSimHoughTransformArray* myHoughTransformArray = (WCSimHoughTransformArray*)(this->HoughTransformArray(myDigitList,myVertex));
00208   
00209   // Find Hough Peak
00210   // ===============
00211   Double_t houghVtxX = myVertex->GetX();
00212   Double_t houghVtxY = myVertex->GetY();
00213   Double_t houghVtxZ = myVertex->GetZ();
00214 
00215   Double_t houghDirX = 0.0;
00216   Double_t houghDirY = 0.0;
00217   Double_t houghDirZ = 0.0;
00218   Double_t houghAngle = 0.0;
00219   Double_t houghHeight = 0.0;
00220  
00221   myHoughTransformArray->FindPeak(houghDirX,houghDirY,houghDirZ,
00222                                   houghAngle,houghHeight);
00223   
00224   // Make New Ring
00225   // =============
00226   std::cout << "  found new ring: (vx,vy,vz)=(" << houghVtxX << "," << houghVtxY << "," << houghVtxZ << ") " << std::endl
00227             << "                  (px,py,pz)=(" << houghDirX << "," << houghDirY << "," << houghDirZ << ") " << std::endl
00228             << "                   angle=" << houghAngle << " height=" << houghHeight << std::endl;
00229 
00230   WCSimRecoRing* myRing = new WCSimRecoRing(houghVtxX,houghVtxY,houghVtxZ,
00231                                             houghDirX,houghDirY,houghDirZ,
00232                                             houghAngle,houghHeight);
00233 
00234   fRingList->push_back(myRing);
00235 
00236   // Return Ring List
00237   // ================
00238   return fRingList;
00239 }
00240 
00241 WCSimHoughTransform* WCSimRingFinder::HoughTransform(WCSimRecoEvent* myEvent, WCSimRecoVertex* myVertex, Double_t myAngle)
00242 {
00243   // get filter digit list
00244   // =====================
00245   std::vector<WCSimRecoDigit*>* myFilterDigitList = (std::vector<WCSimRecoDigit*>*)(myEvent->GetFilterDigitList());
00246 
00247   // run Hough Transform
00248   // ===================
00249   return (WCSimHoughTransform*)(this->HoughTransform(myFilterDigitList,myVertex,myAngle));
00250 }
00251 
00252 WCSimHoughTransform* WCSimRingFinder::HoughTransform(std::vector<WCSimRecoDigit*>* myFilterDigitList, WCSimRecoVertex* myVertex, Double_t myAngle)
00253 {
00254   std::cout << " *** WCSimRingFinder::HoughTransform(...) *** " << std::endl;
00255   std::cout << "  calculate Hough Transform for cone angle: " << myAngle << std::endl;  
00256   
00257   // make new Hough Transform (if necessary)
00258   // =======================================
00259   if( fHoughTransform==0 ){
00260     fHoughTransform = new WCSimHoughTransform(fHoughX,fHoughY);
00261   }
00262 
00263   // reset Hough Transform
00264   // =====================
00265   fHoughTransform->Reset();
00266 
00267   // perform Hough Transform
00268   // =======================
00269   for( UInt_t idigit=0; idigit<myFilterDigitList->size(); idigit++ ){
00270     WCSimRecoDigit* myDigit = (WCSimRecoDigit*)(myFilterDigitList->at(idigit));
00271 
00272     Double_t x = myDigit->GetX();
00273     Double_t y = myDigit->GetY();
00274     Double_t z = myDigit->GetZ();
00275 
00276     Double_t vx = myVertex->GetX();
00277     Double_t vy = myVertex->GetY();
00278     Double_t vz = myVertex->GetZ();
00279 
00280     Double_t omega = 0.0;
00281     Double_t rx = 0.0;
00282     Double_t ry = 0.0;
00283     Double_t rz = 0.0;
00284     Double_t nx = 0.0;
00285     Double_t ny = 0.0;
00286     Double_t nz = 0.0;
00287     Double_t r = 0.0;
00288 
00289     Double_t weight = 1.0;
00290 
00291     for( Int_t n=0; n<fHoughPoints; n++ ){
00292       omega = 360.0*((double)n/(double)fHoughPoints);
00293 
00294       WCSimGeometry::FindCircle(x, y, z,
00295                                 vx, vy, vz,
00296                                 myAngle, omega,
00297                                 rx, ry, rz,
00298                                 nx, ny, nz, r);
00299 
00300       fHoughTransform->Fill(nx,ny,nz,weight);
00301     }
00302   }
00303 
00304   // return result of Hough Transform
00305   // ================================
00306   return fHoughTransform;
00307 }
00308 
00309 WCSimHoughTransformArray* WCSimRingFinder::HoughTransformArray(WCSimRecoEvent* myEvent, WCSimRecoVertex* myVertex)
00310 {
00311   // get filter digit list
00312   // =====================
00313   std::vector<WCSimRecoDigit*>* myFilterDigitList = (std::vector<WCSimRecoDigit*>*)(myEvent->GetFilterDigitList());
00314 
00315   // run Hough Transform
00316   // ===================
00317   return (WCSimHoughTransformArray*)(this->HoughTransformArray(myFilterDigitList,myVertex));
00318 }
00319 
00320 WCSimHoughTransformArray* WCSimRingFinder::HoughTransformArray(std::vector<WCSimRecoDigit*>* myFilterDigitList, WCSimRecoVertex* myVertex)
00321 {
00322   std::cout << " *** WCSimRingFinder::HoughTransformArray(...) *** " << std::endl;
00323   
00324   // make new Hough Transform array (if necessary)
00325   // =============================================
00326   if( fHoughTransformArray==0 ){
00327     fHoughTransformArray = new WCSimHoughTransformArray( fConeAngleBins,fConeAngleMin,fConeAngleMax,
00328                                                          fHoughX,fHoughY);
00329   }
00330 
00331   // reset Hough Transform array
00332   // ===========================
00333   fHoughTransformArray->Reset();
00334 
00335   // perform Hough Transform
00336   // =======================
00337   for( UInt_t idigit=0; idigit<myFilterDigitList->size(); idigit++ ){
00338     WCSimRecoDigit* myDigit = (WCSimRecoDigit*)(myFilterDigitList->at(idigit));
00339 
00340     Double_t x = myDigit->GetX();
00341     Double_t y = myDigit->GetY();
00342     Double_t z = myDigit->GetZ();
00343 
00344     Double_t vx = myVertex->GetX();
00345     Double_t vy = myVertex->GetY();
00346     Double_t vz = myVertex->GetZ();
00347 
00348     Double_t omega = 0.0;
00349     Double_t rx = 0.0;
00350     Double_t ry = 0.0;
00351     Double_t rz = 0.0;
00352     Double_t nx = 0.0;
00353     Double_t ny = 0.0;
00354     Double_t nz = 0.0;
00355     Double_t r = 0.0;
00356 
00357     Double_t weight = 1.0;
00358 
00359     for( Int_t n=0; n<fHoughPoints; n++ ){
00360       omega = 360.0*((double)n/(double)fHoughPoints);
00361 
00362       for( Int_t nAngle=0; nAngle<fHoughTransformArray->GetConeAngleBins(); nAngle++ ){
00363         Double_t myAngle = fHoughTransformArray->GetConeAngle(nAngle);
00364         WCSimHoughTransform* myHoughTransform = (WCSimHoughTransform*)(fHoughTransformArray->GetHoughTransform(nAngle));
00365         
00366         WCSimGeometry::FindCircle(x, y, z,
00367                                   vx, vy, vz,
00368                                   myAngle, omega,
00369                                   rx, ry, rz,
00370                                   nx, ny, nz, r);
00371 
00372         myHoughTransform->Fill(nx,ny,nz,weight);
00373       }
00374     }
00375   }
00376 
00377   // return result of Hough Transform
00378   // ================================
00379   return fHoughTransformArray;
00380 }