00001 #include "DataUtil/PlaneOutline.h"
00002 #include <cmath>
00003 #include <list>
00004 #include <cassert>
00005 #include <algorithm>
00006
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
00034 using Munits::mm;
00035 std::vector< std::pair<float,float> > v;
00036
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
00060 using Munits::mm;
00061 std::vector< std::pair<float,float> > v;
00062
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
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
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
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
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
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
00170 using Munits::inch;
00171 std::vector< std::pair<float,float> > v;
00172
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
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
00188 using Munits::inch;
00189 std::vector< std::pair<float,float> > v;
00190
00191
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
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
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
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
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
00309
00310
00311
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
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
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); }
00350 break;
00351 case PlaneCoverage::kComplete:
00352 result = IsInside(x,y,fFarU);
00353 if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilU); }
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); }
00367 break;
00368 case PlaneCoverage::kComplete:
00369 result = IsInside(x,y,fFarV);
00370 if(result && bCoilOutside){ result = !IsInside(x,y,fFarCoilV); }
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
00385
00386
00387
00388
00389
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
00446
00447
00448
00449
00450
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
00498
00499
00500
00501
00502
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
00543
00544
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
00577
00578
00579
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
00597
00598
00599
00600
00601
00602
00603
00604
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
00619
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
00642
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
00653 for (unsigned int i=0; i<(v.size()-1); i++){
00654
00655
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
00660 }
00661
00662
00663 }
00664 std::list<float>::iterator it = q.begin();
00665 while(it!=q.end()){
00666
00667 it++;
00668 }
00669 }
00670 assert(q.size()%2 == 0);
00671
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;
00690 float u1, v1, u2, v2;
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
00783
00784
00785
00786
00787
00788
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
00808 bool result=false;
00809 float x = ax+bx*s;
00810 float y = ay+by*s;
00811
00812
00813
00814
00815
00816
00817
00818
00819
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
00831
00832
00833
00834
00835
00836
00837
00838
00839
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
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
00957
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
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
01004 htemp=fu;
01005 if(kview==8) htemp=fv;
01006 }
01007 else if(plane<=120){
01008
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
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
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 }