00001
00003
00004
00005
00006
00007
00008
00010
00011 #include "VP1GuideLineSystems/InDetProjHelper.h"
00012 #include "VP1GuideLineSystems/InDetProjParams.h"
00013 #include "VP1Base/VP1Msg.h"
00014
00015
00016
00017
00018 InDetProjHelper * InDetProjHelper::createPixelHelper( IVP1System* system )
00019 {
00020 return new InDetProjHelper( InDetProjParams::surfacethickness(),
00021 InDetProjParams::pixel_data_disttosurface_epsilon(),
00022 InDetProjParams::pixel_barrel_inner_radius(),
00023 InDetProjParams::pixel_barrel_outer_radius(),
00024 InDetProjParams::pixel_barrel_posneg_z(),
00025 InDetProjParams::pixel_endcap_surface_z(),
00026 InDetProjParams::pixel_endcap_surface_length(),
00027 InDetProjParams::pixel_endcap_inner_radius(),
00028 InDetProjParams::pixel_endcap_outer_radius(),
00029 InDetProjParams::pixel_endcap_zasr_innerradius(),
00030 InDetProjParams::pixel_endcap_zasr_endcapz_begin(),
00031 InDetProjParams::pixel_endcap_zasr_squeezefact(),
00032 system );
00033 }
00034
00035
00036 InDetProjHelper * InDetProjHelper::createSCTHelper(IVP1System* system )
00037 {
00038 return new InDetProjHelper( InDetProjParams::surfacethickness(),
00039 InDetProjParams::sct_data_disttosurface_epsilon(),
00040 InDetProjParams::sct_barrel_inner_radius(),
00041 InDetProjParams::sct_barrel_outer_radius(),
00042 InDetProjParams::sct_barrel_posneg_z(),
00043 InDetProjParams::sct_endcap_surface_z(),
00044 InDetProjParams::sct_endcap_surface_length(),
00045 InDetProjParams::sct_endcap_inner_radius(),
00046 InDetProjParams::sct_endcap_outer_radius(),
00047 InDetProjParams::sct_endcap_zasr_innerradius(),
00048 InDetProjParams::sct_endcap_zasr_endcapz_begin(),
00049 InDetProjParams::sct_endcap_zasr_squeezefact(),
00050 system );
00051 }
00052
00053
00054 InDetProjHelper * InDetProjHelper::createTRTHelper(IVP1System* system )
00055 {
00056 return new InDetProjHelper( InDetProjParams::surfacethickness(),
00057 InDetProjParams::trt_data_disttosurface_epsilon(),
00058 InDetProjParams::trt_barrel_inner_radius(),
00059 InDetProjParams::trt_barrel_outer_radius(),
00060 InDetProjParams::trt_barrel_posneg_z(),
00061 InDetProjParams::trt_endcap_surface_z(),
00062 InDetProjParams::trt_endcap_surface_length(),
00063 InDetProjParams::trt_endcap_inner_radius(),
00064 InDetProjParams::trt_endcap_outer_radius(),
00065 InDetProjParams::trt_endcap_zasr_innerradius(),
00066 InDetProjParams::trt_endcap_zasr_endcapz_begin(),
00067 InDetProjParams::trt_endcap_zasr_squeezefact(),
00068 system );
00069 }
00070
00071
00072 class InDetProjHelper::Imp {
00073 public:
00074 InDetProjHelper * theclass;
00075
00076
00077 InDetProjFlags::InDetProjPartsFlags parts;
00078
00079
00080 double surfacethickness;
00081 double data_disttosurface_epsilon;
00082 double barrel_inner_radius;
00083 double barrel_outer_radius;
00084 double barrel_posneg_z;
00085 double endcap_surface_z;
00086 double endcap_surface_length;
00087 double endcap_inner_radius;
00088 double endcap_outer_radius;
00089 double endcap_zasr_innerradius;
00090 double endcap_zasr_endcapz_begin;
00091 double endcap_zasr_squeezefact;
00092
00093
00094
00095 double covercyl_zmin;
00096 double covercyl_zmax;
00097 double covercyl_rmin;
00098 double covercyl_rmax;
00099
00100
00101 void lineCircleIntersection( const HepPoint3D&a, const HepPoint3D&b,const double& r,
00102 double & u1, double& u2 ) const;
00103
00104
00105 void movePoint1ToZPlaneAndPoint2( HepPoint3D& p1, const HepPoint3D& p2, const double& z ) const;
00106 bool clipSegmentToZInterval( HepPoint3D&a, HepPoint3D&b, const double& zmin, const double& zmax ) const;
00107 void movePoint1ToInfiniteCylinderAndPoint2( HepPoint3D&p1, const HepPoint3D&p2, const double& r ) const;
00108 bool clipSegmentToInfiniteHollowCylinder( HepPoint3D&a, HepPoint3D&b,
00109 const double& rmin, const double& rmax,
00110 HepPoint3D&seg2_a, HepPoint3D&seg2_b ) const;
00111
00112 bool clipSegmentToHollowCylinder( HepPoint3D&a, HepPoint3D&b,
00113 const double& rmin, const double& rmax,
00114 const double& zmin, const double& zmax,
00115 HepPoint3D&seg2_a, HepPoint3D&seg2_b ) const;
00116
00117 void clipPathToHollowCylinder( const std::vector<HepPoint3D>& in,
00118 std::set<std::vector<HepPoint3D> >& out,
00119 const double& rmin, const double& rmax,
00120 const double& zmin, const double& zmax ) const;
00121
00122 bool touchesHollowCylinder( const std::vector<HepPoint3D>& path,
00123 const double& rmin, const double& rmax,
00124 const double& zmin, const double& zmax ) const;
00125
00126 void projectPathToInfiniteCylinder( const std::vector<HepPoint3D>& in,
00127 std::set<std::vector<HepPoint3D> >& outset, const double& r ) const;
00128 void projectPathToZPlane( const std::vector<HepPoint3D>& in,
00129 std::set<std::vector<HepPoint3D> >& outset, const double& z ) const;
00130 void projectPathToZPlane_specialZtoR( const std::vector<HepPoint3D>& in,
00131 std::set<std::vector<HepPoint3D> >& outset, const double& z ) const;
00132
00133 };
00134
00135
00136 InDetProjHelper::InDetProjHelper( const double& _surfacethickness,
00137 const double& _data_disttosurface_epsilon,
00138 const double& _barrel_inner_radius,
00139 const double& _barrel_outer_radius,
00140 const double& _barrel_posneg_z,
00141 const double& _endcap_surface_z,
00142 const double& _endcap_surface_length,
00143 const double& _endcap_inner_radius,
00144 const double& _endcap_outer_radius,
00145 const double& _endcap_zasr_innerradius,
00146 const double& _endcap_zasr_endcapz_begin,
00147 const double& _endcap_zasr_squeezefact,
00148 IVP1System* sys )
00149 : VP1HelperClassBase(sys,"InDetProjHelper"), d(new Imp)
00150 {
00151 d->theclass = this;
00152
00153 d->surfacethickness = _surfacethickness;
00154 d->data_disttosurface_epsilon = _data_disttosurface_epsilon;
00155 d->barrel_inner_radius = _barrel_inner_radius;
00156 d->barrel_outer_radius = _barrel_outer_radius;
00157 d->barrel_posneg_z = _barrel_posneg_z;
00158 d->endcap_surface_z = _endcap_surface_z;
00159 d->endcap_surface_length = _endcap_surface_length;
00160 d->endcap_inner_radius = _endcap_inner_radius;
00161 d->endcap_outer_radius = _endcap_outer_radius;
00162 d->endcap_zasr_innerradius = _endcap_zasr_innerradius;
00163 d->endcap_zasr_endcapz_begin = _endcap_zasr_endcapz_begin;
00164 d->endcap_zasr_squeezefact = _endcap_zasr_squeezefact;
00165
00166 d->parts = InDetProjFlags::NoProjections;
00167 d->covercyl_zmin = 0.0;
00168 d->covercyl_zmax = 0.0;
00169 d->covercyl_rmin = 0.0;
00170 d->covercyl_rmax = 0.0;
00171
00172 }
00173
00174
00175 InDetProjHelper::~InDetProjHelper()
00176 {
00177 delete d;
00178 }
00179
00180
00181 InDetProjFlags::InDetProjPartsFlags InDetProjHelper::setParts( InDetProjFlags::InDetProjPartsFlags newparts )
00182 {
00183 if ( d->parts==newparts )
00184 return d->parts;
00185 InDetProjFlags::InDetProjPartsFlags oldparts = d->parts;
00186 d->parts = newparts;
00187
00188
00189 if (d->parts == InDetProjFlags::NoProjections) {
00190 d->covercyl_zmin = 0.0;
00191 d->covercyl_zmax = 0.0;
00192 d->covercyl_rmin = 0.0;
00193 d->covercyl_rmax = 0.0;
00194 return oldparts;
00195 }
00196
00197 bool no_ec_neg = !( d->parts & InDetProjFlags::EndCap_AllNeg );
00198 bool no_ec_pos = !( d->parts & InDetProjFlags::EndCap_AllPos );
00199 bool no_brl_neg = !( d->parts & InDetProjFlags::Barrel_AllNeg );
00200 bool no_brl_pos = !( d->parts & InDetProjFlags::Barrel_AllPos );
00201 bool barrel = d->parts & InDetProjFlags::Barrel_All;
00202 bool endcap = d->parts & InDetProjFlags::EndCap_All;
00203
00204 d->covercyl_zmin = - d->endcap_surface_z - 0.5*d->endcap_surface_length;
00205 if ( no_ec_neg ) {
00206 d->covercyl_zmin = - d->barrel_posneg_z;
00207 if ( no_brl_neg ) {
00208 d->covercyl_zmin = 0.0;
00209 if ( no_brl_pos ) {
00210 d->covercyl_zmin = d->barrel_posneg_z;
00211 if ( no_ec_pos )
00212 d->covercyl_zmin = d->endcap_surface_z + 0.5*d->endcap_surface_length + 1.0e99;
00213 }
00214 }
00215 }
00216 d->covercyl_zmax = d->endcap_surface_z + 0.5*d->endcap_surface_length;
00217 if ( no_ec_pos ) {
00218 d->covercyl_zmax = d->barrel_posneg_z;
00219 if ( no_brl_pos ) {
00220 d->covercyl_zmax = 0.0;
00221 if ( no_brl_neg ) {
00222 d->covercyl_zmax = - d->barrel_posneg_z;
00223 if ( no_ec_neg )
00224 d->covercyl_zmax = - d->endcap_surface_z - 0.5*d->endcap_surface_length - 1.0e99;
00225 }
00226 }
00227 }
00228 if ( d->covercyl_zmin >= d->covercyl_zmax )
00229 d->covercyl_zmin = d->covercyl_zmax = 0;
00230
00231 if ( barrel && endcap ) {
00232 d->covercyl_rmin = std::min(d->barrel_inner_radius,d->endcap_inner_radius);
00233 d->covercyl_rmax = std::max(d->barrel_outer_radius,d->endcap_outer_radius);
00234 } else {
00235 if (barrel) {
00236 d->covercyl_rmin = d->barrel_inner_radius;
00237 d->covercyl_rmax = d->barrel_outer_radius;
00238 } else if (endcap) {
00239 d->covercyl_rmin = d->endcap_inner_radius;
00240 d->covercyl_rmax = d->endcap_outer_radius;
00241 } else {
00242 message("Unforeseen execution path encountered.");
00243 d->covercyl_rmin = 0;
00244 d->covercyl_rmax = 0;
00245 }
00246 }
00247 if ( d->covercyl_rmin >= d->covercyl_rmax )
00248 d->covercyl_rmin = d->covercyl_rmax = 0;
00249 return oldparts;
00250 }
00251
00252
00253 InDetProjFlags::InDetProjPartsFlags InDetProjHelper::parts() const
00254 {
00255 return d->parts;
00256 }
00257
00258
00259 void InDetProjHelper::clipPath( const std::vector<HepPoint3D>& path,
00260 std::set<std::vector<HepPoint3D> >& resulting_subpaths ) const
00261 {
00262 clipPath(path,resulting_subpaths,resulting_subpaths,resulting_subpaths,resulting_subpaths);
00263 }
00264
00265
00266 void InDetProjHelper::clipPath( const std::vector<HepPoint3D>& path,
00267 std::set<std::vector<HepPoint3D> >& resulting_subpaths_barrelA,
00268 std::set<std::vector<HepPoint3D> >& resulting_subpaths_barrelC,
00269 std::set<std::vector<HepPoint3D> >& resulting_subpaths_endcapA,
00270 std::set<std::vector<HepPoint3D> >& resulting_subpaths_endcapC ) const
00271 {
00272 if (verbose())
00273 messageVerbose("clipPath(..) called. Input path has "+QString::number(path.size())+" points.");
00274
00275 resulting_subpaths_barrelA.clear();
00276 resulting_subpaths_barrelC.clear();
00277 resulting_subpaths_endcapA.clear();
00278 resulting_subpaths_endcapC.clear();
00279
00280
00281 if (d->parts == InDetProjFlags::NoProjections ) {
00282 if (verbose())
00283 messageVerbose("All projections currently off.");
00284 return;
00285 }
00286 if ( path.size()<2 ) {
00287 if (verbose())
00288 messageVerbose("Input path too short.");
00289 return;
00290 }
00291
00292
00293
00294
00295
00296 std::set<std::vector<HepPoint3D> > paths_clipped;
00297 d->clipPathToHollowCylinder( path, paths_clipped,
00298 d->covercyl_rmin, d->covercyl_rmax,
00299 d->covercyl_zmin, d->covercyl_zmax );
00300
00301 if (paths_clipped.empty()) {
00302 if (verbose())
00303 messageVerbose("Path entirely outside clip volumes.");
00304 return;
00305 }
00306
00307 const bool enabled_brlA = d->parts & InDetProjFlags::Barrel_AllPos;
00308 const bool enabled_brlC = d->parts & InDetProjFlags::Barrel_AllNeg;
00309 const bool enabled_ecA = d->parts & InDetProjFlags::EndCap_AllPos;
00310 const bool enabled_ecC = d->parts & InDetProjFlags::EndCap_AllNeg;
00311
00312
00313 if ( ( (enabled_brlA?1:0) + (enabled_brlC?1:0) + (enabled_ecA?1:0) + (enabled_ecC?1:0) ) == 1 ) {
00314 if (enabled_brlA) {
00315 resulting_subpaths_barrelA = paths_clipped;
00316 if (verbose())
00317 messageVerbose("clipPath(..) only brlA enabled. Returning.");
00318 return;
00319 }
00320 if (enabled_brlC) {
00321 resulting_subpaths_barrelC = paths_clipped;
00322 if (verbose())
00323 messageVerbose("clipPath(..) only brlC enabled. Returning.");
00324 return;
00325 }
00326 if (enabled_ecA) {
00327 resulting_subpaths_endcapA = paths_clipped;
00328 if (verbose())
00329 messageVerbose("clipPath(..) only ecA enabled. Returning.");
00330 return;
00331 }
00332 if (enabled_ecC) {
00333 resulting_subpaths_endcapC = paths_clipped;
00334 if (verbose())
00335 messageVerbose("clipPath(..) only ecC enabled. Returning.");
00336 return;
00337 }
00338 }
00339
00340
00341
00342
00343
00344 std::set<std::vector<HepPoint3D> >::const_iterator it, itE(paths_clipped.end());
00345 for (it = paths_clipped.begin();it!=itE;++it) {
00346 if ( enabled_brlA )
00347 d->clipPathToHollowCylinder( *it, resulting_subpaths_barrelA, d->barrel_inner_radius, d->barrel_outer_radius, 0, d->barrel_posneg_z );
00348 if ( enabled_brlC )
00349 d->clipPathToHollowCylinder( *it, resulting_subpaths_barrelC, d->barrel_inner_radius, d->barrel_outer_radius, - d->barrel_posneg_z, 0 );
00350 if ( enabled_ecA )
00351 d->clipPathToHollowCylinder( *it, resulting_subpaths_endcapA, d->endcap_inner_radius, d->endcap_outer_radius,
00352 d->endcap_surface_z - d->endcap_surface_length * 0.5, d->endcap_surface_z + d->endcap_surface_length * 0.5 );
00353 if ( enabled_ecC )
00354 d->clipPathToHollowCylinder( *it, resulting_subpaths_endcapC, d->endcap_inner_radius, d->endcap_outer_radius,
00355 - d->endcap_surface_z - d->endcap_surface_length * 0.5, - d->endcap_surface_z + d->endcap_surface_length * 0.5 );
00356 }
00357
00358 messageVerbose("clipPath(..) end.");
00359
00360 }
00361
00362
00363 void InDetProjHelper::Imp::movePoint1ToZPlaneAndPoint2( HepPoint3D& p1, const HepPoint3D& p2, const double& z ) const
00364 {
00365 double dx(p2.x()-p1.x()), dy(p2.y()-p1.y()), dz(p2.z()-p1.z());
00366 if (dz==0.0) {
00367 theclass->message("movePoint1ToZPlaneAndPoint2 Error: Points have same z!!");
00368 return;
00369 }
00370 double s( (z-p1.z())/dz );
00371 p1.set( p1.x()+dx*s, p1.y()+dy*s, z );
00372 }
00373
00374
00375 bool InDetProjHelper::Imp::clipSegmentToZInterval( HepPoint3D&a, HepPoint3D&b,
00376 const double& zmin, const double& zmax ) const
00377 {
00378 if (a.z()<zmin) {
00379 if (b.z()<zmin)
00380 return false;
00381
00382 movePoint1ToZPlaneAndPoint2(a,b,zmin);
00383 if (b.z()>zmax)
00384 movePoint1ToZPlaneAndPoint2(b,a,zmax);
00385 return true;
00386 } else {
00387 if (b.z()<zmin) {
00388
00389 movePoint1ToZPlaneAndPoint2(b,a,zmin);
00390 if (a.z()>zmax)
00391 movePoint1ToZPlaneAndPoint2(a,b,zmax);
00392 return true;
00393 } else {
00394
00395 if (a.z()>zmax) {
00396 if (b.z()>zmax)
00397 return false;
00398 movePoint1ToZPlaneAndPoint2(a,b,zmax);
00399 return true;
00400 } else {
00401
00402 if (b.z()>zmax)
00403 movePoint1ToZPlaneAndPoint2(b,a,zmax);
00404 return true;
00405 }
00406 }
00407 }
00408 }
00409
00410
00411 void InDetProjHelper::Imp::movePoint1ToInfiniteCylinderAndPoint2( HepPoint3D&p1, const HepPoint3D&p2,const double& r ) const
00412 {
00413
00414
00415
00416
00417 double p1r(p1.r());
00418 double dr(p2.r()-p1r);
00419 if (dr==0.0) {
00420 theclass->message("movePoint1ToInfiniteCylinderAndPoint2 Error: Points have same r!!");
00421 return;
00422 }
00423 double s((r-p1r)/dr);
00424 double t(1.0-s);
00425 p1.set( p1.x()*t + p2.x()*s, p1.y()*t + p2.y()*s, p1.z()*t + p2.z()*s );
00426
00427 }
00428
00429
00430 void InDetProjHelper::Imp::lineCircleIntersection( const HepPoint3D&a, const HepPoint3D&b,const double& r,
00431 double & u1, double& u2 ) const
00432 {
00433 const double dx = b.x()-a.x();
00434 const double dy = b.y()-a.y();
00435 double A = dx*dx+dy*dy;
00436 if (A==0.0) {
00437
00438 u1 = u2 = ( a.x()*a.x()+a.y()*a.y() == r*r ? 0.0 : 1.0e99 );
00439 return;
00440 }
00441 double B = 2.0*( a.x()*dx + a.y()*dy );
00442 double C = a.x()*a.x()+a.y()*a.y() - r*r;
00443 double D = B*B-4*A*C;
00444
00445 if (D>0.0) {
00446
00447 double sqrtD = sqrt(D);
00448 u1 = 0.5 * ( -B - sqrtD) / A;
00449 u2 = 0.5 * ( -B + sqrtD) / A;
00450 } else if (D<0.0) {
00451
00452 u1 = u2 = -1.0e99;
00453 } else {
00454
00455 u1 = u2 = -0.5*B/A;
00456 }
00457
00458 }
00459
00460
00461 bool InDetProjHelper::Imp::clipSegmentToInfiniteHollowCylinder( HepPoint3D&a, HepPoint3D&b,
00462 const double& rmin, const double& rmax,
00463 HepPoint3D&seg2_a, HepPoint3D&seg2_b ) const
00464 {
00465
00466
00467
00468
00469
00470
00471
00472
00473
00474
00475
00476
00477
00478
00479
00480 const double ar2 = a.x()*a.x()+a.y()*a.y();
00481 const double br2 = b.x()*b.x()+b.y()*b.y();
00482 const double rmin2 = rmin*rmin;
00483
00484 if (ar2 <= rmin2 && br2 <= rmin2 ) {
00485
00486
00487
00488 return false;
00489 }
00490
00491 if ( (a.x()<=-rmax&&b.x()<=-rmax) || (a.x()>=rmax&&b.x()>=rmax) || (a.y()<=-rmax&&b.y()<=-rmax)|| (a.y()>=rmax&&b.y()>=rmax) ) {
00492
00493
00494
00495 return false;
00496 }
00497
00498
00499 const double dx = b.x()-a.x();
00500 const double dy = b.y()-a.y();
00501 const double rmax2 = rmax*rmax;
00502 if (dx==0.0&&dy==0.0) {
00503
00504
00505
00506 return ar2<=rmax2;
00507 }
00508
00509 const double u = - (a.y()*dy+a.x()*dx)/(dx*dx+dy*dy);
00510 const double px = ( u <= 0 ? a.x() : ( u >= 1 ? b.x() : a.x()+u*dx ) );
00511 const double py = ( u <= 0 ? a.y() : ( u >= 1 ? b.y() : a.y()+u*dy ) );
00512 const double pr2 = px*px+py*py;
00513 if (pr2>=rmax2) {
00514
00515
00516
00517 return false;
00518
00519
00520 }
00521
00522 seg2_a=seg2_b;
00523
00524 if (pr2>=rmin2&&ar2<=rmax2&&br2<=rmax2) {
00525
00526
00527
00528
00529
00530 return true;
00531 }
00532
00533
00534 if (ar2>rmax2||br2>rmax2) {
00535
00536
00537 double u1, u2;
00538 lineCircleIntersection(a,b,rmax,u1,u2);
00539 if (u1==u2) {
00540
00541 theclass->message("This should never happen(1).");
00542
00543 return false;
00544 }
00545 HepPoint3D asave(a);
00546 if (u1>0&&u1<1) {
00547
00548 a = a+u1*(b-a);
00549
00550
00551 }
00552 if (u2>0&&u2<1) {
00553
00554 b = asave+u2*(b-asave);
00555
00556
00557 }
00558 }
00559
00560 if (pr2>=rmin2) {
00561
00562
00563 return true;
00564 }
00565
00566 double u1, u2;
00567 lineCircleIntersection(a,b,rmin,u1,u2);
00568
00569 if (u1>0&&u1<1) {
00570 if (u2>0&&u2<1) {
00571
00572
00573
00574 seg2_b = b;
00575 b = a+u1*(seg2_b-a);
00576 seg2_a=a+u2*(seg2_b-a);
00577
00578
00579 return true;
00580 }
00581 b = a+u1*(b-a);
00582
00583 return true;
00584 }
00585 if (u2>0&&u2<1)
00586 a = a+u2*(b-a);
00587
00588 return true;
00589 }
00590
00591
00592
00593
00594 bool InDetProjHelper::Imp::clipSegmentToHollowCylinder( HepPoint3D&a, HepPoint3D&b,
00595 const double& rmin, const double& rmax,
00596 const double& zmin, const double& zmax,
00597 HepPoint3D&seg2_a, HepPoint3D&seg2_b ) const
00598 {
00599
00600
00601
00602
00603
00604
00605
00606
00607
00608
00609 if (!clipSegmentToZInterval(a,b,zmin,zmax)) {
00610
00611
00612
00613 return false;
00614 }
00615
00616
00617
00618
00619
00620
00621
00622
00623
00624 if (!clipSegmentToInfiniteHollowCylinder(a,b,rmin,rmax,seg2_a,seg2_b)) {
00625
00626
00627
00628 return false;
00629 }
00630
00631
00632
00633
00634
00635
00636
00637
00638
00639
00640
00641
00642
00643 return true;
00644 }
00645
00646
00647 void InDetProjHelper::Imp::clipPathToHollowCylinder( const std::vector<HepPoint3D>& in,
00648 std::set<std::vector<HepPoint3D> >& out,
00649 const double& rmin, const double& rmax,
00650 const double& zmin, const double& zmax ) const
00651 {
00652
00653
00654
00655
00656
00657
00658
00659
00660 out.clear();
00661 if (rmin>=rmax||rmin<0||zmin>=zmax) {
00662 theclass->message("clipPathToHollowCylinder Error: Non-sensical cylinder parameters!");
00663 return;
00664 }
00665 const unsigned n=in.size();
00666 if (n<2)
00667 return;
00668
00669 HepPoint3D a,b;
00670 HepPoint3D seg2_a,seg2_b;
00671 std::vector<HepPoint3D> v;
00672 for (unsigned i = 1; i<n; ++i) {
00673
00674
00675 a = in.at(i-1);
00676 b = in.at(i);
00677 if ( clipSegmentToHollowCylinder( a,b,rmin,rmax,zmin,zmax,seg2_a,seg2_b ) ) {
00678 if (v.empty()) {
00679 v.push_back(a);
00680 v.push_back(b);
00681 if (seg2_a!=seg2_b) {
00682 out.insert(v);
00683 v.clear();
00684 v.push_back(seg2_a);
00685 v.push_back(seg2_b);
00686 }
00687 } else {
00688
00689
00690 if ( v.back() != a ) {
00691 theclass->messageDebug("ERROR: Inconsistency found while building clip part");
00692 out.insert(v);
00693 v.clear();
00694 v.push_back(a);
00695 }
00696 v.push_back(b);
00697 if (seg2_a!=seg2_b) {
00698 out.insert(v);
00699 v.clear();
00700 v.push_back(seg2_a);
00701 v.push_back(seg2_b);
00702 }
00703 }
00704 } else {
00705
00706
00707 if (!v.empty()) {
00708 out.insert(v);
00709 v.clear();
00710 }
00711 }
00712 }
00713 if (!v.empty())
00714 out.insert(v);
00715 }
00716
00717
00718 void InDetProjHelper::Imp::projectPathToInfiniteCylinder( const std::vector<HepPoint3D>& in,
00719 std::set<std::vector<HepPoint3D> >& outset, const double& r ) const
00720 {
00721 std::vector<HepPoint3D> out(in);
00722 std::vector<HepPoint3D>::iterator it(out.begin()), itE(out.end());
00723 double s;
00724 for (;it!=itE;++it) {
00725 if ( it->x()==0.0 && it->y()==0.0 ) {
00726 theclass->message("ERROR: Point has x==0 and y==0. Ambiguous projection of point.");
00727 it->setX(1.0);
00728 }
00729 s = r / sqrt( it->x()*it->x()+it->y()*it->y() );
00730 it->setX(it->x()*s);
00731 it->setY(it->y()*s);
00732 }
00733 outset.insert(out);
00734 }
00735
00736
00737 void InDetProjHelper::Imp::projectPathToZPlane( const std::vector<HepPoint3D>& in,
00738 std::set<std::vector<HepPoint3D> >& outset, const double& z ) const
00739 {
00740 std::vector<HepPoint3D> out(in);
00741 std::vector<HepPoint3D>::iterator it(out.begin()), itE(out.end());
00742 for (;it!=itE;++it)
00743 it->setZ(z);
00744 outset.insert(out);
00745 }
00746
00747
00748
00749 void InDetProjHelper::transformECPointToZPlane_specialZtoR( HepPoint3D& p,
00750 const double& planeZ,
00751 const double& planeRBegin,
00752 const double& endcapZBegin,
00753 const double& squeezeFactor )
00754 {
00755 if ( p.x()==0.0 && p.y()==0.0 ) {
00756 VP1Msg::message("InDetProjHelper::transformECPointToZPlane_specialZtoR ERROR: "
00757 "Point has x==0 and y==0. Ambiguous projection of point.");
00758 p.setX(1.0);
00759 }
00760 const double r = planeRBegin + (fabs(p.z())-endcapZBegin)/squeezeFactor;
00761 const double s = r / sqrt( p.x()*p.x()+p.y()*p.y() );
00762 p.setX(p.x()*s);
00763 p.setY(p.y()*s);
00764 p.setZ(planeZ);
00765 }
00766
00767
00768 void InDetProjHelper::Imp::projectPathToZPlane_specialZtoR( const std::vector<HepPoint3D>& in,
00769 std::set<std::vector<HepPoint3D> >& outset,
00770 const double& z ) const
00771 {
00772 std::vector<HepPoint3D> out(in);
00773 std::vector<HepPoint3D>::iterator it(out.begin()), itE(out.end());
00774 for (;it!=itE;++it)
00775 InDetProjHelper::transformECPointToZPlane_specialZtoR(*it,
00776 z,
00777 endcap_zasr_innerradius,
00778 endcap_zasr_endcapz_begin,
00779 endcap_zasr_squeezefact);
00780 outset.insert(out);
00781 }
00782
00783
00784 void InDetProjHelper::projectPath( const std::vector<HepPoint3D>& path,
00785 std::set<std::vector<HepPoint3D> >& resulting_projs ) const
00786 {
00787 projectPath(path,resulting_projs,resulting_projs,resulting_projs,resulting_projs);
00788 }
00789
00790
00791 void InDetProjHelper::projectPath( const std::vector<HepPoint3D>& path,
00792 std::set<std::vector<HepPoint3D> >& resulting_projections_barrelA,
00793 std::set<std::vector<HepPoint3D> >& resulting_projections_barrelC,
00794 std::set<std::vector<HepPoint3D> >& resulting_projections_endcapA,
00795 std::set<std::vector<HepPoint3D> >& resulting_projections_endcapC ) const
00796 {
00797 if (verbose())
00798 messageVerbose("projectPath(..) called. Input path has "+QString::number(path.size())+" points.");
00799
00800 resulting_projections_barrelA.clear();
00801 resulting_projections_barrelC.clear();
00802 resulting_projections_endcapA.clear();
00803 resulting_projections_endcapC.clear();
00804
00805
00806 if (d->parts == InDetProjFlags::NoProjections ) {
00807 if (verbose())
00808 messageVerbose("All projections currently off.");
00809 return;
00810 }
00811 if ( path.size()<2 ) {
00812 if (verbose())
00813 messageVerbose("Input path too short.");
00814 return;
00815 }
00816
00817
00818
00819 std::set<std::vector<HepPoint3D> > paths_brlA, paths_brlC, paths_ecA,paths_ecC;
00820 clipPath( path,paths_brlA, paths_brlC, paths_ecA,paths_ecC);
00821
00822
00823
00824
00825
00826 const double eps = d->data_disttosurface_epsilon;
00827 const double endcapeps(-5*mm);
00828
00829 std::set<std::vector<HepPoint3D> >::const_iterator it,itE;
00830
00831 if (d->parts & InDetProjFlags::Barrel_AllPos) {
00832 itE = paths_brlA.end();
00833 if ( d->parts & InDetProjFlags::BarrelCentral )
00834 for ( it = paths_brlA.begin(); it!=itE; ++it )
00835 d->projectPathToZPlane( *it, resulting_projections_barrelA, 0.5*d->surfacethickness+eps );
00836 if ( d->parts & InDetProjFlags::BarrelPositive )
00837 for ( it = paths_brlA.begin(); it!=itE; ++it )
00838 d->projectPathToZPlane( *it, resulting_projections_barrelA, d->barrel_posneg_z - eps );
00839 }
00840 if ( d->parts & InDetProjFlags::Barrel_AllNeg ) {
00841 itE = paths_brlC.end();
00842 if ( d->parts & InDetProjFlags::BarrelCentral )
00843 for ( it = paths_brlC.begin(); it!=itE; ++it )
00844 d->projectPathToZPlane( *it, resulting_projections_barrelC, - 0.5*d->surfacethickness - eps);
00845 if ( d->parts & InDetProjFlags::BarrelNegative )
00846 for ( it = paths_brlC.begin(); it!=itE; ++it )
00847 d->projectPathToZPlane( *it, resulting_projections_barrelC, - d->barrel_posneg_z );
00848 }
00849 if ( d->parts & InDetProjFlags::EndCap_AllPos ) {
00850 itE = paths_ecA.end();
00851 if ( d->parts & InDetProjFlags::EndCapInnerPositive )
00852 for ( it = paths_ecA.begin(); it!=itE; ++it )
00853 d->projectPathToInfiniteCylinder( *it, resulting_projections_endcapA, d->endcap_inner_radius + eps+endcapeps );
00854 if ( d->parts & InDetProjFlags::EndCapOuterPositive )
00855 for ( it = paths_ecA.begin(); it!=itE; ++it )
00856 d->projectPathToInfiniteCylinder( *it, resulting_projections_endcapA, d->endcap_outer_radius + eps+endcapeps );
00857
00858 if ( d->parts & InDetProjFlags::TRT_EndCapZToRCentral )
00859 for ( it = paths_ecA.begin(); it!=itE; ++it )
00860 d->projectPathToZPlane_specialZtoR( *it, resulting_projections_endcapA,
00861 0.5*d->surfacethickness + eps );
00862
00863 if ( d->parts & InDetProjFlags::TRT_EndCapZToRPositive )
00864 for ( it = paths_ecA.begin(); it!=itE; ++it )
00865 d->projectPathToZPlane_specialZtoR( *it, resulting_projections_endcapA,
00866 d->barrel_posneg_z - 0.5*d->surfacethickness - eps );
00867 }
00868 if ( d->parts & InDetProjFlags::EndCap_AllNeg ) {
00869 itE = paths_ecC.end();
00870 if ( d->parts & InDetProjFlags::EndCapInnerNegative )
00871 for ( it = paths_ecC.begin(); it!=itE; ++it )
00872 d->projectPathToInfiniteCylinder( *it, resulting_projections_endcapC, d->endcap_inner_radius + eps+endcapeps );
00873 if ( d->parts & InDetProjFlags::EndCapOuterNegative )
00874 for ( it = paths_ecC.begin(); it!=itE; ++it )
00875 d->projectPathToInfiniteCylinder( *it, resulting_projections_endcapC, d->endcap_outer_radius + eps+endcapeps );
00876
00877 if ( d->parts & InDetProjFlags::TRT_EndCapZToRCentral )
00878 for ( it = paths_ecC.begin(); it!=itE; ++it )
00879 d->projectPathToZPlane_specialZtoR( *it, resulting_projections_endcapC,
00880 - 0.5*d->surfacethickness - eps );
00881
00882 if ( d->parts & InDetProjFlags::TRT_EndCapZToRNegative )
00883 for ( it = paths_ecC.begin(); it!=itE; ++it )
00884 d->projectPathToZPlane_specialZtoR( *it, resulting_projections_endcapC,
00885 - d->barrel_posneg_z + 0.5*d->surfacethickness + eps );
00886 }
00887
00888 }
00889
00890
00891 InDetProjHelper::PartsFlags InDetProjHelper::touchedParts( const std::vector<HepPoint3D>& path ) const
00892 {
00893 if (verbose())
00894 messageVerbose("touchedParts(..) called. Input path has "+QString::number(path.size())+" points.");
00895 PartsFlags touchedparts = NoParts;
00896 if ( d->touchesHollowCylinder(path,d->barrel_inner_radius, d->barrel_outer_radius, 0, d->barrel_posneg_z) )
00897 touchedparts |= BarrelA;
00898 if ( d->touchesHollowCylinder(path,d->barrel_inner_radius, d->barrel_outer_radius, - d->barrel_posneg_z, 0) )
00899 touchedparts |= BarrelC;
00900 if ( d->touchesHollowCylinder(path,d->endcap_inner_radius, d->endcap_outer_radius,
00901 d->endcap_surface_z - d->endcap_surface_length * 0.5, d->endcap_surface_z + d->endcap_surface_length * 0.5 ) )
00902 touchedparts |= EndCapA;
00903 if ( d->touchesHollowCylinder(path, d->endcap_inner_radius, d->endcap_outer_radius,
00904 - d->endcap_surface_z - d->endcap_surface_length * 0.5, - d->endcap_surface_z + d->endcap_surface_length * 0.5) )
00905 touchedparts |= EndCapC;
00906 return touchedparts;
00907 }
00908
00909
00910 bool InDetProjHelper::Imp::touchesHollowCylinder( const std::vector<HepPoint3D>& path,
00911 const double& rmin, const double& rmax,
00912 const double& zmin, const double& zmax ) const
00913 {
00914 const double rmin2(rmin*rmin), rmax2(rmax*rmax);
00915 double r2;
00916 std::vector<HepPoint3D>::const_iterator it(path.begin()), itE(path.end());
00917 for (;it!=itE;++it) {
00918 if (it->z()<zmin)
00919 continue;
00920 if (it->z()>zmax)
00921 continue;
00922 r2 = it->x()*it->x()+it->y()*it->y();
00923 if (r2<rmin2)
00924 continue;
00925 if (r2<=rmax2)
00926 return true;
00927 }
00928 return false;
00929 }