loon: PlaneOutline.cxx Source File

00001 #include "DataUtil/PlaneOutline.h"
00002 #include <cmath>
00003 #include <list>
00004 #include <cassert>
00005 #include <algorithm>
00006 //#include <iostream>
00007 
00008 #include "TPolyLine.h"
00009 #include "TMath.h"
00010 #include "TH2.h"
00011 #include "TExec.h"
00012 #include "TList.h"
00013 
00014 #include "Conventions/Munits.h"
00015 
00016 std::vector< std::pair<float, float> > PlaneOutline::fNearPartialV  = PlaneOutline::InitNearPartialV();
00017 std::vector< std::pair<float, float> > PlaneOutline::fNearPartialU  = PlaneOutline::InitNearPartialU();
00018 std::vector< std::pair<float, float> > PlaneOutline::fNearFullV     = PlaneOutline::InitNearFullV();
00019 std::vector< std::pair<float, float> > PlaneOutline::fNearFullU     = PlaneOutline::InitNearFullU();
00020 std::vector< std::pair<float, float> > PlaneOutline::fNearFullCoilV = PlaneOutline::InitNearFullCoilV();
00021 std::vector< std::pair<float, float> > PlaneOutline::fNearFullCoilU = PlaneOutline::InitNearFullCoilU();
00022 std::vector< std::pair<float, float> > PlaneOutline::fNearSteel     = PlaneOutline::InitNearSteel();
00023 
00024 std::vector< std::pair<float, float> > PlaneOutline::fFarV          = PlaneOutline::InitFarV();
00025 std::vector< std::pair<float, float> > PlaneOutline::fFarU          = PlaneOutline::InitFarU();
00026 std::vector< std::pair<float, float> > PlaneOutline::fFarCoilV      = PlaneOutline::InitFarCoilV();
00027 std::vector< std::pair<float, float> > PlaneOutline::fFarCoilU      = PlaneOutline::InitFarCoilU();
00028 std::vector< std::pair<float, float> > PlaneOutline::fFarSteel      = PlaneOutline::InitFarSteel();
00029 
00030 PlaneOutline::PlaneOutline() {}
00031 
00032 std::vector< std::pair<float,float> > PlaneOutline::InitFarCoilU(){
00033   // coil hole of full u plane
00034   using Munits::mm;
00035   std::vector< std::pair<float,float> > v;
00036   // input in u,v coordinates
00037   v.push_back( std::pair<float,float>(41*mm, -300*mm) );
00038   v.push_back( std::pair<float,float>(243*mm, -92*mm) );
00039   v.push_back( std::pair<float,float>(243*mm, 92*mm) );
00040   v.push_back( std::pair<float,float>(41*mm, 300*mm) );
00041   v.push_back( std::pair<float,float>(-41*mm, 300*mm) );
00042   v.push_back( std::pair<float,float>(-243*mm, 92*mm) );
00043   v.push_back( std::pair<float,float>(-243*mm, -92*mm) );
00044   v.push_back( std::pair<float,float>(-41*mm, -300*mm) );
00045   v.push_back( std::pair<float,float>(41*mm, -300*mm) );
00046   
00047   for (unsigned int i=0; i<v.size(); i++){    
00048     float x=sqrt(1.0/2.0)*(v[i].first - v[i].second);
00049     float y=sqrt(1.0/2.0)*(v[i].first + v[i].second);
00050     v[i].first = x;
00051     v[i].second = y;
00052   }
00053   
00054 
00055   return v;
00056 }
00057 
00058 std::vector< std::pair<float,float> > PlaneOutline::InitFarCoilV(){
00059   // coil hole of full u plane
00060   using Munits::mm;
00061   std::vector< std::pair<float,float> > v;
00062   // input in u,v coordinates
00063 
00064   v.push_back( std::pair<float,float>(-300*mm, 41*mm) );
00065   v.push_back( std::pair<float,float>(-92*mm, 243*mm) );
00066   v.push_back( std::pair<float,float>(92*mm, 243*mm) );
00067   v.push_back( std::pair<float,float>(300*mm, 41*mm) );
00068   v.push_back( std::pair<float,float>(300*mm, -41*mm) );
00069   v.push_back( std::pair<float,float>(92*mm, -243*mm) );
00070   v.push_back( std::pair<float,float>(-92*mm, -243*mm) );
00071   v.push_back( std::pair<float,float>(-300*mm, -41*mm) );
00072   v.push_back( std::pair<float,float>(-300*mm, 41*mm) );
00073   
00074   for (unsigned int i=0; i<v.size(); i++){    
00075     float x=sqrt(1.0/2.0)*(v[i].first - v[i].second);
00076     float y=sqrt(1.0/2.0)*(v[i].first + v[i].second);
00077     v[i].first = x;
00078     v[i].second = y;
00079   }
00080   
00081 
00082   return v;
00083 }
00084 
00085 
00086 std::vector< std::pair<float,float> > PlaneOutline::InitNearFullCoilU(){
00087   // coil hole of full u plane
00088   using Munits::inch;
00089   std::vector< std::pair<float,float> > v;
00090   v.push_back( std::pair<float,float>(-12.22*inch, 12.56*inch) );
00091   v.push_back( std::pair<float,float>(2.89*inch, 15.10*inch) );
00092   v.push_back( std::pair<float,float>(14.93*inch, 3.06*inch) );
00093   v.push_back( std::pair<float,float>(12.22*inch, -12.22*inch) );
00094   v.push_back( std::pair<float,float>(-3.06*inch, -14.93*inch) );
00095   v.push_back( std::pair<float,float>(-14.93*inch, -3.22*inch) );
00096   v.push_back( std::pair<float,float>(-12.22*inch, 12.56*inch) );
00097 
00098   return v;
00099 }
00100 
00101 
00102 std::vector< std::pair<float,float> > PlaneOutline::InitNearFullCoilV(){
00103   // coil hole of full v plane
00104   using Munits::inch;
00105   std::vector< std::pair<float,float> > v;
00106   v.push_back( std::pair<float,float>(-3.13*inch, 14.60*inch) );
00107   v.push_back( std::pair<float,float>(12.16*inch, 12.34*inch) );
00108   v.push_back( std::pair<float,float>(14.25*inch, -3.30*inch) );
00109   v.push_back( std::pair<float,float>(1.74*inch, -14.77*inch) );
00110   v.push_back( std::pair<float,float>(-12.69*inch, -12.34*inch) );
00111   v.push_back( std::pair<float,float>(-14.77*inch, 2.78*inch) );
00112   v.push_back( std::pair<float,float>(-3.13*inch, 14.60*inch) );
00113 
00114   return v;
00115 }
00116 
00117 
00118 std::vector< std::pair<float,float> > PlaneOutline::InitNearPartialV(){
00119   // partial v plane
00120   using Munits::inch;
00121   std::vector< std::pair<float,float> > v;
00122   v.push_back( std::pair<float,float>(-1.96*inch, 12.72*inch) );
00123   v.push_back( std::pair<float,float>(58.00*inch, 73.15*inch) );
00124   v.push_back( std::pair<float,float>(103.75*inch, 27.89*inch) );
00125   v.push_back( std::pair<float,float>(108.15*inch, 32.04*inch) );
00126   v.push_back( std::pair<float,float>(108.15*inch, -32.07*inch) );
00127   v.push_back( std::pair<float,float>(62.90*inch, -77.58*inch) );
00128   v.push_back( std::pair<float,float>(-1.79*inch, -12.50*inch) );
00129   v.push_back( std::pair<float,float>(10.05*inch, 0.0*inch) );
00130   v.push_back( std::pair<float,float>(-1.96*inch, 12.72*inch) );
00131 
00132   return v;
00133 }
00134 
00135 std::vector< std::pair<float,float> > PlaneOutline::InitNearPartialU(){
00136   // partial u plane
00137   using Munits::inch;
00138   std::vector< std::pair<float,float> > v;
00139   v.push_back( std::pair<float,float>(-1.79*inch, 12.30*inch) );
00140   v.push_back( std::pair<float,float>(62.92*inch, 77.46*inch) );
00141   v.push_back( std::pair<float,float>(108.36*inch, 32.00*inch) );
00142   v.push_back( std::pair<float,float>(108.36*inch, -32.03*inch) );
00143   v.push_back( std::pair<float,float>(103.88*inch, -27.55*inch) );
00144   v.push_back( std::pair<float,float>(58.45*inch, -72.98*inch) );
00145   v.push_back( std::pair<float,float>(-1.79*inch, -12.54*inch) );
00146   v.push_back( std::pair<float,float>(10.09*inch, 0.0*inch) );
00147   v.push_back( std::pair<float,float>(-1.79*inch, 12.30*inch) );
00148 
00149   return v;
00150 }
00151 
00152 std::vector< std::pair<float,float> > PlaneOutline::InitNearFullV(){
00153   // full v plane
00154   using Munits::inch;
00155   std::vector< std::pair<float,float> > v;
00156   v.push_back( std::pair<float,float>(-67.52*inch, 4.93*inch) );
00157   v.push_back( std::pair<float,float>(4.17*inch, 76.10*inch) );
00158   v.push_back( std::pair<float,float>(72.18*inch, 72.18*inch) );
00159   v.push_back( std::pair<float,float>(113.11*inch, 31.52*inch) );
00160   v.push_back( std::pair<float,float>(113.11*inch, -31.80*inch) );
00161   v.push_back( std::pair<float,float>(71.66*inch, -73.77*inch) );
00162   v.push_back( std::pair<float,float>(9.63*inch, -72.45*inch) );
00163   v.push_back( std::pair<float,float>(-67.52*inch, 4.93*inch) );
00164 
00165   return v;
00166 }
00167 
00168 std::vector< std::pair<float,float> > PlaneOutline::InitNearFullU(){
00169   // full u plane
00170   using Munits::inch;
00171   std::vector< std::pair<float,float> > v;
00172   //  v.push_back( std::pair<float,float>(-68.22*inch, 5.36*inch) );
00173   v.push_back( std::pair<float,float>(-68.22*inch, -5.36*inch) );
00174   v.push_back( std::pair<float,float>(9.67*inch, 73.44*inch) );
00175   v.push_back( std::pair<float,float>(72.02*inch, 73.44*inch) );
00176   v.push_back( std::pair<float,float>(114.28*inch, 32.08*inch) );
00177   v.push_back( std::pair<float,float>(114.28*inch, -32.08*inch) );
00178   v.push_back( std::pair<float,float>(73.04*inch, -73.59*inch) );
00179   v.push_back( std::pair<float,float>(4.07*inch, -76.88*inch) );
00180   //  v.push_back( std::pair<float,float>(-68.22*inch, 5.36*inch) );
00181   v.push_back( std::pair<float,float>(-68.22*inch, -5.36*inch) );
00182 
00183   return v;
00184 }
00185 
00186 std::vector< std::pair<float,float> > PlaneOutline::InitNearSteel(){
00187   // near steel plane
00188   using Munits::inch;
00189   std::vector< std::pair<float,float> > v;
00190 
00191   // Steel pln (coords about the center, *not* MINOS coords)
00192   v.push_back( std::pair<float,float>(-121.43*inch,   0.47*inch) );
00193   v.push_back( std::pair<float,float>(-121.43*inch,   9.28*inch) );
00194   v.push_back( std::pair<float,float>( -69.80*inch,  60.91*inch) );
00195   v.push_back( std::pair<float,float>( -69.80*inch,  75.04*inch) );
00196   v.push_back( std::pair<float,float>(  69.80*inch,  75.04*inch) );
00197   v.push_back( std::pair<float,float>(  69.80*inch,  60.91*inch) );
00198   v.push_back( std::pair<float,float>( 121.43*inch,   9.28*inch) );
00199   v.push_back( std::pair<float,float>( 121.43*inch,   0.47*inch) );
00200   v.push_back( std::pair<float,float>( 120.25*inch,   0.47*inch) );
00201   v.push_back( std::pair<float,float>( 117.30*inch,   0.47*inch) );
00202   v.push_back( std::pair<float,float>( 110.16*inch,   3.42*inch) );
00203   v.push_back( std::pair<float,float>(  95.24*inch,  -1.48*inch) );
00204   v.push_back( std::pair<float,float>(  95.24*inch, -35.47*inch) );
00205   v.push_back( std::pair<float,float>(  69.80*inch, -60.91*inch) );
00206   v.push_back( std::pair<float,float>(  69.80*inch, -75.04*inch) );
00207   v.push_back( std::pair<float,float>( -69.80*inch, -75.04*inch) );
00208   v.push_back( std::pair<float,float>( -69.80*inch, -60.91*inch) );
00209   v.push_back( std::pair<float,float>( -95.24*inch, -35.47*inch) );
00210   v.push_back( std::pair<float,float>( -95.24*inch,  -1.48*inch) );
00211   v.push_back( std::pair<float,float>(-110.24*inch,   3.42*inch) );
00212   v.push_back( std::pair<float,float>(-115.14*inch,   3.42*inch) );
00213   v.push_back( std::pair<float,float>(-120.25*inch,   0.47*inch) );
00214   v.push_back( std::pair<float,float>(-121.43*inch,   0.47*inch) );
00215 
00216   // Correct to MINOS coords
00217   const double xOffset = +0.5578;
00218   for (unsigned int i = 0; i < v.size(); i++){
00219     v[i].first += xOffset;
00220   }
00221 
00222   return v;
00223 }
00224 
00225 std::vector< std::pair<float,float> > PlaneOutline::InitFarSteel(){
00226   // far steel plane
00227   using Munits::cm;
00228   std::vector< std::pair<float,float> > v;
00229 
00230   v.push_back( std::pair<float,float>(  143.7724*cm, -421.0036*cm) );
00231   v.push_back( std::pair<float,float>(  139.4268*cm, -421.0036*cm) );
00232   v.push_back( std::pair<float,float>(  116.8712*cm, -398.8619*cm) );
00233   v.push_back( std::pair<float,float>( -116.9622*cm, -398.8619*cm) );
00234   v.push_back( std::pair<float,float>( -139.1040*cm, -421.0036*cm) );
00235   v.push_back( std::pair<float,float>( -143.4496*cm, -421.0036*cm) );
00236   v.push_back( std::pair<float,float>( -421.5665*cm, -142.8867*cm) );
00237   v.push_back( std::pair<float,float>( -421.5665*cm, -137.7134*cm) );
00238   v.push_back( std::pair<float,float>( -399.8386*cm, -115.9855*cm) );
00239   v.push_back( std::pair<float,float>( -399.8386*cm,   84.7388*cm) );
00240   v.push_back( std::pair<float,float>( -406.6674*cm,   94.2577*cm) );
00241   v.push_back( std::pair<float,float>( -429.8438*cm,  101.9141*cm) );
00242   v.push_back( std::pair<float,float>( -439.9834*cm,  101.9141*cm) );
00243   v.push_back( std::pair<float,float>( -457.3658*cm,   94.2577*cm) );
00244   v.push_back( std::pair<float,float>( -457.3658*cm,   94.2577*cm) );
00245   v.push_back( std::pair<float,float>( -457.3658*cm,  109.3637*cm) );
00246   v.push_back( std::pair<float,float>( -143.4496*cm,  423.2799*cm) );
00247   v.push_back( std::pair<float,float>( -139.1040*cm,  423.2799*cm) );
00248   v.push_back( std::pair<float,float>( -116.9622*cm,  401.1382*cm) );
00249   v.push_back( std::pair<float,float>(  117.2850*cm,  401.1382*cm) );
00250   v.push_back( std::pair<float,float>(  139.4268*cm,  423.2799*cm) );
00251   v.push_back( std::pair<float,float>(  143.7724*cm,  423.2799*cm) );
00252   v.push_back( std::pair<float,float>(  457.6886*cm,  109.3637*cm) );
00253   v.push_back( std::pair<float,float>(  457.6886*cm,   94.2577*cm) );
00254   v.push_back( std::pair<float,float>(  454.7915*cm,   94.2577*cm) );
00255   v.push_back( std::pair<float,float>(  447.1350*cm,  101.5003*cm) );
00256   v.push_back( std::pair<float,float>(  429.1319*cm,  101.5003*cm) );
00257   v.push_back( std::pair<float,float>(  406.9902*cm,   94.0507*cm) );
00258   v.push_back( std::pair<float,float>(  400.7823*cm,   87.6358*cm) );
00259   v.push_back( std::pair<float,float>(  400.1614*cm, -115.9855*cm) );
00260   v.push_back( std::pair<float,float>(  421.8893*cm, -138.3342*cm) );
00261   v.push_back( std::pair<float,float>(  421.8893*cm, -142.8876*cm) );
00262   v.push_back( std::pair<float,float>(  143.7724*cm, -421.0036*cm) );
00263 
00264   return v;
00265 } 
00266 
00267 std::vector< std::pair<float,float> > PlaneOutline::InitFarU(){
00268   // full u plane
00269   const float d=4.0*Munits::meter;
00270   const float L=2.0*d*tan(TMath::Pi()/8.0);
00271   std::vector< std::pair<float,float> > v;
00272   v.push_back( std::pair<float,float>(-L/2.0, d) );
00273   v.push_back( std::pair<float,float>(L/2.0, d) );
00274   v.push_back( std::pair<float,float>(d, L/2.0) );
00275   v.push_back( std::pair<float,float>(d, -L/2.0) );
00276   v.push_back( std::pair<float,float>(L/2.0, -d) );
00277   v.push_back( std::pair<float,float>(-L/2.0, -d) );
00278   v.push_back( std::pair<float,float>(-d, -L/2.0) );
00279   v.push_back( std::pair<float,float>(-d, L/2.0) );
00280   v.push_back( std::pair<float,float>(-L/2.0, d) );
00281 
00282   return v;
00283 }
00284 
00285 std::vector< std::pair<float,float> > PlaneOutline::InitFarV(){
00286   // full u plane
00287   const float d=4.0*Munits::meter;
00288   const float L=2.0*d*tan(TMath::Pi()/8.0);
00289   std::vector< std::pair<float,float> > v;
00290   v.push_back( std::pair<float,float>(-L/2.0, d) );
00291   v.push_back( std::pair<float,float>(L/2.0, d) );
00292   v.push_back( std::pair<float,float>(d, L/2.0) );
00293   v.push_back( std::pair<float,float>(d, -L/2.0) );
00294   v.push_back( std::pair<float,float>(L/2.0, -d) );
00295   v.push_back( std::pair<float,float>(-L/2.0, -d) );
00296   v.push_back( std::pair<float,float>(-d, -L/2.0) );
00297   v.push_back( std::pair<float,float>(-d, L/2.0) );
00298   v.push_back( std::pair<float,float>(-L/2.0, d) );
00299 
00300   return v;
00301 }
00302 
00303 
00304 
00305 bool PlaneOutline::IsInside(float xp, float yp, 
00306                             const std::vector< std::pair<float,float> >& v) const
00307 {
00308   // Is (xp,yp) inside the plane outline represented by v?
00309 
00310 
00311   // taken from TMath::IsInside and modified to use my data structure
00312   float xint;
00313   int inter = 0;
00314   int np=v.size(); 
00315    for (int i=0;i<np-1;i++) {
00316       if (v[i].second == v[i+1].second) continue;
00317       if (yp <= v[i].second && yp <= v[i+1].second) continue;
00318       if (v[i].second < yp && v[i+1].second < yp) continue;
00319       xint = v[i].first 
00320         + (yp-v[i].second)*(v[i+1].first-v[i].first)/(v[i+1].second-v[i].second);
00321       if (xp < xint) inter++;
00322    }
00323    if (inter%2) return true;
00324    return false;
00325 }
00326 
00327 bool PlaneOutline::IsInsideSteel(float x, float y, Detector::Detector_t det) const {
00328   // Is the (x,y) coord inside the outline of a steel plane in detector det (near/far)?
00329 
00330   bool result = false;
00331   if (det == Detector::kNear) result = IsInside(x,y,fNearSteel);
00332   if (det == Detector::kFar)  result = IsInside(x,y,fFarSteel);
00333 
00334   return result;
00335 }
00336 
00337 bool PlaneOutline::IsInside(float x, float y, PlaneView::PlaneView_t pv, 
00338                             PlaneCoverage::PlaneCoverage_t pc, bool bCoilOutside) const {
00339   // Is 'x,y' inside the outline of the plane type represented by 'pv' and 'pc'?
00340 
00341   bool result=false;
00342   if( pv==PlaneView::kU ){
00343     switch(pc){
00344     case PlaneCoverage::kNearPartial:
00345       result = IsInside(x,y,fNearPartialU);
00346       break;
00347     case PlaneCoverage::kNearFull:
00348       result = IsInside(x,y,fNearFullU);
00349       if(result && bCoilOutside){ result = !IsInside(x,y,fNearFullCoilU); } // coil hole region
00350       break;
00351     case PlaneCoverage::kComplete:
00352       result = IsInside(x,y,fFarU);
00353       if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilU); } // coil hole region      
00354       break;
00355     default:
00356       break;
00357     }
00358   }
00359   else if (pv==PlaneView::kV){
00360     switch(pc){
00361     case PlaneCoverage::kNearPartial:
00362       result = IsInside(x,y,fNearPartialV);
00363       break;
00364     case PlaneCoverage::kNearFull:
00365       result = IsInside(x,y,fNearFullV);
00366       if(result && bCoilOutside){ result = !IsInside(x,y,fNearFullCoilV); } // coil hole region
00367       break;
00368     case PlaneCoverage::kComplete:
00369       result = IsInside(x,y,fFarV);
00370       if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilV); } // coil hole region      
00371       break;
00372     default:
00373       break;      
00374     }
00375   }
00376   return result;
00377  
00378 }
00379 
00380 
00381 bool PlaneOutline::DistanceToEdge(float x, float y, PlaneView::PlaneView_t pv,
00382                                   PlaneCoverage::PlaneCoverage_t pc, 
00383                                   float& dist, float& xedge, float& yedge) const {
00384   // compute the distance from 'x,y' to the closest point on the edge 
00385   // of the plane denoted by the 'pc' and 'pv' arguments
00386   //
00387   // if a distance can be computed then set 'dist' to the minimum distance
00388   // and 'xedge', 'yedge' to the closest point on that edge and return true
00389   // otherwise return false
00390 
00391  
00392   bool result=false;
00393   if( pv==PlaneView::kU ){
00394     float dist2, x2,y2;
00395     switch(pc){
00396     case PlaneCoverage::kNearPartial:
00397       dist = DistanceToEdge(x,y,fNearPartialU,xedge,yedge);
00398       result=true;
00399       break;
00400     case PlaneCoverage::kNearFull:
00401       dist = DistanceToEdge(x,y,fNearFullU,xedge,yedge);
00402       dist2 = DistanceToEdge(x,y,fNearFullCoilU,x2,y2);
00403       if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
00404       result=true;
00405       break;
00406     case PlaneCoverage::kComplete:
00407       dist = DistanceToEdge(x,y,fFarU,xedge,yedge);
00408       dist2 = DistanceToEdge(x,y,fFarCoilU,x2,y2);
00409       if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
00410       result=true;
00411       break;
00412     default:
00413       break;
00414     }
00415   }
00416   else if (pv==PlaneView::kV){
00417     float dist2, x2, y2;
00418     switch(pc){
00419     case PlaneCoverage::kNearPartial:
00420       dist = DistanceToEdge(x,y,fNearPartialV,xedge,yedge);
00421       result=true;
00422       break;
00423     case PlaneCoverage::kNearFull:
00424       dist = DistanceToEdge(x,y,fNearFullV,xedge,yedge);
00425       dist2 = DistanceToEdge(x,y,fNearFullCoilV,x2,y2);
00426       if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
00427       result=true;
00428       break;
00429     case PlaneCoverage::kComplete:
00430       dist = DistanceToEdge(x,y,fFarV,xedge,yedge);
00431       dist2 = DistanceToEdge(x,y,fFarCoilV,x2,y2);
00432       if(dist2<dist) {yedge=y2; xedge=x2; dist=dist2;}
00433       result=true;
00434       break;
00435     default:
00436       break;      
00437     }
00438   }
00439   return result; 
00440 }
00441 
00442 bool PlaneOutline::DistanceToOuterEdge(float x, float y, PlaneView::PlaneView_t pv,
00443                                   PlaneCoverage::PlaneCoverage_t pc, 
00444                                   float& dist, float& xedge, float& yedge) const {
00445   // compute the distance from 'x,y' to the closest point on the edge 
00446   // of the plane denoted by the 'pc' and 'pv' arguments
00447   //
00448   // if a distance can be computed then set 'dist' to the minimum distance
00449   // and 'xedge', 'yedge' to the closest point on that edge and return true
00450   // otherwise return false
00451 
00452  
00453   bool result=false;
00454   if( pv==PlaneView::kU ){
00455     switch(pc){
00456     case PlaneCoverage::kNearPartial:
00457       dist = DistanceToEdge(x,y,fNearPartialU,xedge,yedge);
00458       result=true;
00459       break;
00460     case PlaneCoverage::kNearFull:
00461       dist = DistanceToEdge(x,y,fNearFullU,xedge,yedge);
00462       result=true;
00463       break;
00464     case PlaneCoverage::kComplete:
00465       dist = DistanceToEdge(x,y,fFarU,xedge,yedge);
00466       result=true;
00467       break;
00468     default:
00469       break;
00470     }
00471   }
00472   else if (pv==PlaneView::kV){
00473     switch(pc){
00474     case PlaneCoverage::kNearPartial:
00475       dist = DistanceToEdge(x,y,fNearPartialV,xedge,yedge);
00476       result=true;
00477       break;
00478     case PlaneCoverage::kNearFull:
00479       dist = DistanceToEdge(x,y,fNearFullV,xedge,yedge);
00480       result=true;
00481       break;
00482     case PlaneCoverage::kComplete:
00483       dist = DistanceToEdge(x,y,fFarV,xedge,yedge);
00484       result=true;
00485       break;
00486     default:
00487       break;      
00488     }
00489   }
00490   return result; 
00491 }
00492 
00493 
00494 bool PlaneOutline::DistanceToInnerEdge(float x, float y, PlaneView::PlaneView_t pv,
00495                                   PlaneCoverage::PlaneCoverage_t pc, 
00496                                   float& dist, float& xedge, float& yedge) const {
00497   // compute the distance from 'x,y' to the closest point on the edge 
00498   // of the plane denoted by the 'pc' and 'pv' arguments
00499   //
00500   // if a distance can be computed then set 'dist' to the minimum distance
00501   // and 'xedge', 'yedge' to the closest point on that edge and return true
00502   // otherwise return false
00503 
00504  
00505   bool result=false;
00506   if( pv==PlaneView::kU ){
00507     switch(pc){
00508     case PlaneCoverage::kNearFull:
00509       dist = DistanceToEdge(x,y,fNearFullCoilU,xedge,yedge);
00510       result=true;
00511       break;
00512     case PlaneCoverage::kComplete:
00513       dist = DistanceToEdge(x,y,fFarCoilU,xedge,yedge);
00514       result=true;
00515       break;
00516     default:
00517       break;
00518     }
00519   }
00520   else if (pv==PlaneView::kV){
00521     switch(pc){
00522     case PlaneCoverage::kNearFull:
00523       dist = DistanceToEdge(x,y,fNearFullCoilV,xedge,yedge);
00524       result=true;
00525       break;
00526     case PlaneCoverage::kComplete:
00527       dist = DistanceToEdge(x,y,fFarCoilV,xedge,yedge);
00528       result=true;
00529       break;
00530     default:
00531       break;      
00532     }
00533   }
00534   return result; 
00535 }
00536 
00537 
00538 float PlaneOutline::DistanceToEdge(float x, float y, 
00539                                    const std::vector< std::pair<float, float> >& v,
00540                                    float& xedge, float& yedge) const{
00541 
00542   // compute the distance from x,y to the closest point on the edge 
00543   // of the outline represented in v
00544   // works both inside and outside the 
00545   float save_dist = 0., save_x = 0., save_y = 0.;
00546   for (unsigned int i=1; i<v.size(); i++){
00547     float dist = ClosestPointToSegment(v[i-1].first, v[i-1].second,
00548                                     v[i].first, v[i].second,
00549                                     x,y, xedge, yedge);
00550     if(i==1){
00551       save_dist=dist;
00552       save_x=xedge;
00553       save_y=yedge;
00554     }
00555     else if(dist<save_dist){
00556       save_dist=dist;
00557       save_x=xedge;
00558       save_y=yedge;
00559     }
00560   }
00561   xedge=save_x;
00562   yedge=save_y;
00563   return save_dist;
00564 
00565 }
00566 
00567 
00568 float PlaneOutline::ClosestPointToLine(const float& x1, 
00569                                        const float& y1,
00570                                        const float& x2,
00571                                        const float& y2,
00572                                        const float& xpt,
00573                                        const float& ypt,
00574                                        float& x, float& y){
00575   
00576   // given an infinite line defined by points (x1,y1) (x2,y2) 
00577   // and a spatial point (xpt,ypt)
00578   // report the closest point on the line (x,y) to the spatial point
00579   // and compute the distance from (xpt,ypt) to (x,y)
00580   float u= (xpt-x1)*(x2-x1)+(ypt-y1)*(y2-y1);
00581   u/=sqrt( pow(x2-x1,2) + pow(y2-y1,2));
00582   x=x1+u*(x2-x1);
00583   y=y1+u*(y2-y1);
00584   float dist =sqrt(pow(xpt-x, 2)+pow(ypt-y,2));
00585   return dist;
00586 }
00587 
00588 
00589 float PlaneOutline::ClosestPointToSegment(const float& x1, 
00590                                           const float& y1,
00591                                           const float& x2,
00592                                           const float& y2,
00593                                           const float& xpt,
00594                                           const float& ypt,
00595                                           float& x, float& y){
00596   // given a line segment with end points (x1,y1) (x2,y2) 
00597   // and a spatial point (xpt,ypt)
00598   // report the closest point on the line (x,y) to the spatial point
00599   // and compute the distance from (xpt,ypt) to (x,y)
00600   //
00601   // if V(s) is a vector pointing to the point a distance s along the line
00602   // and W is the vector pointing to (xpt,ypt)
00603   // then (x,y) can be found by varying s to minimize |V-W|
00604   // i.e. d|V(s)-W|/ds = 0
00605   //
00606   
00607 
00608   float dist=0.0;
00609 
00610   float d = sqrt(pow(x2-x1,2)+pow(y2-y1,2) );
00611   if(d==0){
00612     x=x1;y=y1;
00613   }
00614   else {
00615     float nx = (x2-x1)/d;  float ny = (y2-y1)/d;
00616     float s = ( (xpt - x1)*nx + (ypt - y1)*ny )/(nx*nx + ny*ny);
00617     
00618     // check that s refers to a point inside the line segment
00619     // if not, set x,y to an end point
00620     if(s<0){
00621       x=x1;
00622       y=y1;
00623     }
00624     else if(s>d){
00625       x=x2;
00626       y=y2;
00627     }
00628     else{
00629       x=x1+s*nx;
00630       y=y1+s*ny;
00631     }
00632 
00633   }
00634   dist=sqrt( pow(x-xpt,2)+pow(y-ypt,2));
00635   return dist;
00636 }
00637 
00638 bool PlaneOutline::DepthInView(float x1, float y1, float x2, float y2, float& d,
00639                                std::vector< std::pair<float,float> >& v) const{
00640   d=0.0;
00641   // find all intersections
00642   // store their distance along segment defined by (x1,y1),(x2,y2)
00643   std::list<float> q;
00644   for (unsigned int i=0; i<(v.size()-1); i++){
00645     float s,l;
00646     if(WhereIntersect(x1,y1,x2,y2, v[i].first, v[i].second, 
00647                       v[i+1].first, v[i+1].second,s,l) ) {
00648       q.push_back(s);
00649     }
00650   }
00651   if(q.size()%2 !=0 ){
00652     //    std::cerr<<"["<<x1<<", "<<y1<<"]   ["<<x2<<", "<<y2<<"]"<<std::endl;
00653     for (unsigned int i=0; i<(v.size()-1); i++){
00654       //      std::cerr<<"("<<v[i].first<<", "<<v[i].second<<") --> ("
00655       //               <<v[i+1].first<<", "<<v[i+1].second<<")";
00656       float s,l;
00657       if(WhereIntersect(x1,y1,x2,y2, v[i].first, v[i].second, 
00658                         v[i+1].first, v[i+1].second,s,l) ) {
00659         //      std::cerr<<" : (s,l) = ("<<s<<", "<<l<<")";
00660       }
00661       //      std::cerr<<std::endl;
00662 
00663     }
00664     std::list<float>::iterator it = q.begin();
00665     while(it!=q.end()){
00666       //      std::cerr<<(*it)<<std::endl;
00667       it++;
00668     }    
00669   }
00670   assert(q.size()%2 == 0); // must be an even number of intersections
00671   //  std::sort(q.begin(),q.end());
00672   q.sort();
00673   std::list<float>::iterator it = q.begin();
00674   while(it!=q.end()){
00675     float a=*it;
00676     it++;
00677     float b=*it;
00678     assert( (b-a) > 0);
00679     d += (b-a);
00680     it++;
00681   }
00682   return true;
00683 }
00684 
00685 float PlaneOutline::DepthInView(float d, PlaneView::PlaneView_t myview,
00686                                 PlaneView::PlaneView_t pv, PlaneCoverage::PlaneCoverage_t pc) const {
00687 
00688   const float is2 = 1.0/sqrt(2.0);
00689   float x1, y1, x2, y2; // will be used later
00690   float u1, v1, u2, v2; // for bookkeeping, defined here due to scoping
00691   switch(myview){
00692   case PlaneView::kX:
00693     y1=d; y2=d;
00694     x1=-8.2*Munits::meter; x2=8.2*Munits::meter;
00695     break;
00696   case PlaneView::kY:
00697     x1=d; x2=d;
00698     y1=-8.2*Munits::meter; y2=8.2*Munits::meter;
00699     break;
00700   case PlaneView::kU:
00701     v1=d; v2=d;
00702     u1=-8.2*Munits::meter; u2=8.2*Munits::meter;
00703     x1 = is2*(u1-v1);
00704     x2 = is2*(u2-v2);
00705     y1 = is2*(u1+v1);
00706     y2 = is2*(u2+v2);
00707     break;
00708   case PlaneView::kV:
00709     u1=d; u2=d;
00710     v1=-8.2*Munits::meter; v2=8.2*Munits::meter;
00711     x1 = is2*(u1-v1);
00712     x2 = is2*(u2-v2);
00713     y1 = is2*(u1+v1);
00714     y2 = is2*(u2+v2);
00715     break;
00716   default:
00717     return 0;
00718     break;
00719   }
00720   
00721   bool result=false;
00722   float dist=0.0;
00723   float dist2=0.0;
00724   if( pv==PlaneView::kU ){
00725     switch(pc){
00726     case PlaneCoverage::kNearPartial:
00727       if(DepthInView(x1,y1,x2,y2,dist,fNearPartialU)) result=true;
00728       break;
00729     case PlaneCoverage::kNearFull:
00730       if(DepthInView(x1,y1,x2,y2,dist,fNearFullU)) {
00731         if(DepthInView(x1,y1,x2,y2,dist2, fNearFullCoilU) ){
00732           dist-=dist2;
00733           result=true;
00734         }
00735       }
00736       break;
00737     case PlaneCoverage::kComplete:
00738       if(DepthInView(x1,y1,x2,y2,dist,fFarU)) {
00739         if(DepthInView(x1,y1,x2,y2,dist2, fFarCoilU) ){
00740           dist-=dist2;
00741           result=true;
00742         }
00743       }
00744       break;
00745     default:
00746       break;
00747     }
00748   }
00749   else if( pv==PlaneView::kV ){
00750     switch(pc){
00751     case PlaneCoverage::kNearPartial:
00752       if(DepthInView(x1,y1,x2,y2,dist,fNearPartialV)) result=true;
00753       break;
00754     case PlaneCoverage::kNearFull:
00755       if(DepthInView(x1,y1,x2,y2,dist,fNearFullV)) {
00756         if(DepthInView(x1,y1,x2,y2,dist2, fNearFullCoilV) ){
00757           dist-=dist2;
00758           result=true;
00759         }
00760       }
00761       break;
00762     case PlaneCoverage::kComplete:
00763       if(DepthInView(x1,y1,x2,y2,dist,fFarV)) {
00764         if(DepthInView(x1,y1,x2,y2,dist2, fFarCoilV) ){
00765           dist-=dist2;
00766           result=true;
00767         }
00768       }
00769       break;
00770     default:
00771       break;
00772     }
00773   }
00774   return dist; 
00775 }
00776 
00777 
00778 bool PlaneOutline::WhereIntersect(float x1, float y1, float x2, float y2,
00779                                   float x3, float y3, float x4, float y4,
00780                                   float& s, float& l){
00781 
00782   // x = ax + s*bx   =   cx + l*dx
00783   // y = ay + s*by   =   cy + l*dy
00784 
00785   // solve for s and l the distance along the two infinite lines, defined by
00786   // (x1,y1),(x2,y2) and (x3,y3),(x4,y4), where the lines intersect
00787   // 
00788   // if these two line *segments* intersect then return true
00789 
00790   float b=sqrt( pow(x2-x1,2) + pow(y2-y1,2) );
00791   float bx = (x2-x1)/b;  float by = (y2-y1)/b;
00792 
00793   float d=sqrt( pow(x4-x3,2) + pow(y4-y3,2) );
00794   float dx = (x4-x3)/d;  float dy = (y4-y3)/d;
00795 
00796   float ax = x1; float ay = y1;
00797   float cx = x3; float cy = y3;
00798 
00799   float den=dx*by-bx*dy;
00800   if(den!=0){
00801     s = (dx*(cy-ay) + dy*(ax-cx))/den;
00802     if(fabs(dx)>0) l = ((ax - cx) + s*bx)/dx;
00803     else  l = ((ay - cy) + s*by)/dy;
00804   }
00805   else return false;
00806 
00807   // is the intersection point within the line segments?
00808   bool result=false;
00809   float x = ax+bx*s;
00810   float y = ay+by*s;
00811   //  std::cout<<"(x,y)=("<<x<<", "<<y<<")"<<std::endl;
00812 
00813   /*
00814   if( ( (x>=x1 && x<=x2) || (x>=x2 && x<=x1) ) &&
00815       ( (y>=y1 && y<=y2) || (y>=y2 && y<=y1) ) &&
00816       ( (x>=x3 && x<=x4) || (x>=x4 && x<=x3) ) &&
00817       ( (y>=y3 && y<=y4) || (y>=y4 && y<=y3) ) ){
00818 
00819     result=true;
00820   }
00821   */
00822 
00823   bool t1 = ( (x>=x1 && x<=x2) || (x>=x2 && x<=x1) || (fabs(x-x2)<1e-6 && fabs(x-x1)<1e-6 && fabs(x1-x2)<1e-6));
00824   bool t2 = ( (y>=y1 && y<=y2) || (y>=y2 && y<=y1) || (fabs(y-y2)<1e-6 && fabs(y-y1)<1e-6 && fabs(y1-y2)<1e-6));
00825   bool t3 = ( (x>=x3 && x<=x4) || (x>=x4 && x<=x3) || (fabs(x-x4)<1e-6 && fabs(x-x3)<1e-6 && fabs(x3-x4)<1e-6));
00826   bool t4 = ( (y>=y3 && y<=y4) || (y>=y4 && y<=y3) || (fabs(y-y4)<1e-6 && fabs(y-y3)<1e-6 && fabs(y3-y4)<1e-6));
00827   
00828 
00829   /*
00830   std::cout<<"[";
00831   if(t1)std::cout<<"yes ";
00832   else std::cout<<"no ";
00833   if(t2)std::cout<<"yes ";
00834   else std::cout<<"no ";
00835   if(t3)std::cout<<"yes ";
00836   else std::cout<<"no ";
00837   if(t4)std::cout<<"yes ";
00838   else std::cout<<"no ";
00839   std::cout<<"]"<<std::endl;
00840   */
00841   if(t1 && t2 && t3 && t4) result=true;
00842 
00843   return result;
00844 }
00845 
00846 
00847 void PlaneOutline::GetSteelOutline(Detector::Detector_t det,
00848                                    TPolyLine*& outer       ){
00849   std::vector< std::pair<float,float> >* v = (det == Detector::kNear)? &fNearSteel : &fFarSteel;
00850 
00851   outer = new TPolyLine(v->size());
00852   for (unsigned int i=0; i<v->size(); i++){
00853     outer->SetPoint(i,(*v)[i].first,(*v)[i].second);
00854   }
00855 
00856   return;
00857 }
00858 
00859 void PlaneOutline::GetOutline(PlaneView::PlaneView_t pv,
00860                               PlaneCoverage::PlaneCoverage_t pc, 
00861                               TPolyLine*& outer, TPolyLine*& inner){
00862   
00863   std::vector< std::pair<float,float> >* v = &fNearFullU;
00864   std::vector< std::pair<float,float> >* vin = &fNearFullCoilU;
00865   bool result=true;
00866   bool doin=false;
00867   if( pv==PlaneView::kU ){
00868     switch(pc){
00869     case PlaneCoverage::kNearPartial:
00870       v=&fNearPartialU;
00871       break;
00872     case PlaneCoverage::kNearFull:
00873       //      std::cout<<"hello"<<std::endl;
00874       v=&fNearFullU;
00875       vin=&fNearFullCoilU; doin=true;
00876       break;
00877     case PlaneCoverage::kComplete:
00878       v=&fFarU;
00879       vin=&fFarCoilU; doin=true;
00880       break;
00881     default:
00882       result=false;
00883       break;
00884     }
00885   }
00886   else if (pv==PlaneView::kV){
00887     switch(pc){
00888     case PlaneCoverage::kNearPartial:
00889       v=&fNearPartialV;
00890       break;
00891     case PlaneCoverage::kNearFull:
00892       v=&fNearFullV;
00893       vin=&fNearFullCoilV; doin=true;
00894       break;
00895     case PlaneCoverage::kComplete:
00896       v=&fFarV;
00897       vin=&fFarCoilV; doin=true;
00898       break;
00899     default:
00900       result=false;
00901       break;    
00902     }
00903   }
00904   
00905   outer=0;
00906   inner=0;
00907   if(result){
00908     outer = new TPolyLine(v->size());
00909     for (unsigned int i=0; i<v->size(); i++){
00910       outer->SetPoint(i,(*v)[i].first,(*v)[i].second);
00911     }
00912   }
00913   if(doin){
00914     inner = new TPolyLine(vin->size());
00915     for (unsigned int i=0; i<vin->size(); i++){
00916       inner->SetPoint(i,(*vin)[i].first,(*vin)[i].second);
00917     }
00918   }
00919 
00920 }
00921 
00922 
00923 TH2* PlaneOutline::DepthInViewHist(PlaneCoverage::PlaneCoverage_t pc, 
00924                                    PlaneView::PlaneView_t pv){
00925   
00926   TH2F* h = new TH2F("h", ";distance (m); view direction;", 
00927                      410, -4.1*Munits::meter, 4.1*Munits::meter, 
00928                      9, -0.5, 8.5);
00929   h->GetYaxis()->SetBinLabel(1, ""); 
00930   h->GetYaxis()->SetBinLabel(2, "X"); 
00931   h->GetYaxis()->SetBinLabel(3, ""); 
00932   h->GetYaxis()->SetBinLabel(4, "Y"); 
00933   h->GetYaxis()->SetBinLabel(5, ""); 
00934   h->GetYaxis()->SetBinLabel(6, "U"); 
00935   h->GetYaxis()->SetBinLabel(7, ""); 
00936   h->GetYaxis()->SetBinLabel(8, "V"); 
00937   h->GetYaxis()->SetBinLabel(9, ""); 
00938   
00939   
00940  for(int i=1; i<=h->GetNbinsX(); i++){
00941    float s = h->GetXaxis()->GetBinCenter(i);
00942    float depth = DepthInView(s,PlaneView::kX,pv,pc);     
00943    h->Fill(s,1,depth);
00944    depth = DepthInView(s,PlaneView::kY,pv,pc); 
00945    h->Fill(s,3,depth);
00946    depth = DepthInView(s,PlaneView::kU,pv,pc);     
00947    h->Fill(s,5,depth);
00948    depth = DepthInView(s,PlaneView::kV,pv,pc);     
00949    h->Fill(s,7,depth);   
00950  }
00951  h->SetDirectory(0);
00952  return h;
00953 }
00954 
00955 TH2* PlaneOutline::GetNDPlanesHist(PlaneView::PlaneView_t view){
00956   // view here is the direction you are looking
00957   // so, to see V readout view=U
00958   if(!(view==PlaneView::kU || view==PlaneView::kV)) return 0;
00959 
00960   TH2* pu = DepthInViewHist(PlaneCoverage::kNearPartial,
00961                             PlaneView::kU);
00962   
00963   TH2* fu = DepthInViewHist(PlaneCoverage::kNearFull,
00964                             PlaneView::kU);
00965   
00966   TH2* pv = DepthInViewHist(PlaneCoverage::kNearPartial,
00967                             PlaneView::kV);
00968   
00969   TH2* fv = DepthInViewHist(PlaneCoverage::kNearFull,
00970                             PlaneView::kV);
00971   
00972   double offset=0.000;
00973   double lowlim=offset;
00974   int n=282;
00975   double pitch=0.0594;
00976   double highlim=n*pitch;
00977 
00978   
00979   TH2* h =0;
00980   if(view==PlaneView::kU) {
00981     h = new TH2F("po_vvz","Transverse vs Z view - V Planes",
00982                  n,lowlim,highlim,
00983                  fv->GetNbinsX(), 
00984                  fv->GetXaxis()->GetBinLowEdge(1),
00985                  fv->GetXaxis()->GetBinUpEdge(fv->GetNbinsX()));
00986   }
00987   else if(view==PlaneView::kV) {
00988     h = new TH2F("po_uvz","Transverse vs Z view - U Planes",
00989                  n,lowlim,highlim,
00990                  fu->GetNbinsX(), 
00991                  fu->GetXaxis()->GetBinLowEdge(1),
00992                  fu->GetXaxis()->GetBinUpEdge(fu->GetNbinsX()));
00993   }
00994 
00995   // kview: 6=U, 8=V
00996   int kview=6;
00997   int pstart=2;
00998   if(view==PlaneView::kV) { kview=8; pstart=1;}
00999   
01000   for(int plane=pstart; plane<n; plane+=2){
01001     TH1* htemp=0;
01002     if((plane-1)%5==0){
01003       // full
01004       htemp=fu;
01005       if(kview==8) htemp=fv;
01006     }
01007     else if(plane<=120){
01008       // partial
01009       htemp=pu;
01010       if(kview==8) htemp=pv;
01011     }
01012     if(htemp!=0){
01013       for(int j=1; j<=htemp->GetNbinsX(); j++){
01014         
01015         float bw=htemp->GetBinContent(j,kview);
01016         float tpos=htemp->GetXaxis()->GetBinCenter(j);
01017         int bin = h->GetYaxis()->FindBin(tpos);
01018         if(bw>0) bw=0.5;
01019         h->SetBinContent(plane+1,bin,bw);
01020       }
01021     }
01022   }
01023   
01024   int coil_bin_high=h->GetYaxis()->FindBin(0.5);
01025   int coil_bin_low=h->GetYaxis()->FindBin(-0.5);
01026   
01027   for(int plane=2; plane<n; plane+=1){
01028     //  for(int j=coil_bin_low; j<=coil_bin_high; j++){
01029     h->SetBinContent(plane+1,coil_bin_low,0.95);
01030     h->SetBinContent(plane+1,coil_bin_high,0.95);
01031     //  }
01032     
01033   }
01034   
01035   h->SetMaximum(1.0);
01036   h->SetMinimum(0.0);
01037   //h->Draw("col");
01038 
01039   TExec* ex = 0;
01040   if(view==6)
01041     ex = new TExec("ex_vvz","gStyle->SetPalette(8,0);");
01042   else
01043     ex = new TExec("ex_uvz","gStyle->SetPalette(8,0);");
01044   
01045   h->GetListOfFunctions()->Add(ex);
01046   
01047   delete pu;
01048   delete fu;
01049   delete pv;
01050   delete fv;
01051   
01052   h->SetDirectory(0);
01053 
01054   return h;
01055 }