CMS 3D CMS Logo

BladeShapeBuilderFromDet.cc
Go to the documentation of this file.
2 
4 
7 
8 #include <iomanip>
9 
10 using namespace std;
11 
12 
13 namespace {
14 
15  pair<DiskSectorBounds *, GlobalVector>
16  computeBounds( const vector<const GeomDet*>& dets,
17  const Plane& plane)
18  {
19  Surface::PositionType tmpPos = dets.front()->surface().position();
20 
21 
22  float rmin(plane.toLocal(tmpPos).perp());
23  float rmax(plane.toLocal(tmpPos).perp());
24  float zmin(plane.toLocal(tmpPos).z());
25  float zmax(plane.toLocal(tmpPos).z());
26  float phimin(plane.toLocal(tmpPos).phi());
27  float phimax(plane.toLocal(tmpPos).phi());
28 
29  for(vector<const GeomDet*>::const_iterator it=dets.begin(); it!=dets.end(); it++){
30  vector<GlobalPoint> corners = BoundingBox().corners( (*it)->specificSurface() );
31 
32  for(vector<GlobalPoint>::const_iterator i=corners.begin();
33  i!=corners.end(); i++) {
34 
35  float r = plane.toLocal(*i).perp();
36  float z = plane.toLocal(*i).z();
37  float phi = plane.toLocal(*i).phi();
38  rmin = min( rmin, r);
39  rmax = max( rmax, r);
40  zmin = min( zmin, z);
41  zmax = max( zmax, z);
42  if ( PhiLess()( phi, phimin)) phimin = phi;
43  if ( PhiLess()( phimax, phi)) phimax = phi;
44  }
45  // in addition to the corners we have to check the middle of the
46  // det +/- length/2, since the min (max) radius for typical fw
47  // dets is reached there
48 
49  float rdet = (*it)->position().perp();
50  float height = (*it)->surface().bounds().width();
51  rmin = min( rmin, rdet-height/2.F);
52  rmax = max( rmax, rdet+height/2.F);
53 
54 
55  }
56 
57  if (!PhiLess()(phimin, phimax))
58  edm::LogError("TkDetLayers") << " BladeShapeBuilderFromDet : "
59  << "Something went wrong with Phi Sorting !" ;
60  float zPos = (zmax+zmin)/2.;
61  float phiWin = phimax - phimin;
62  float phiPos = (phimax+phimin)/2.;
63  float rmed = (rmin+rmax)/2.;
64  if ( phiWin < 0. ) {
65  if ( (phimin < Geom::pi() / 2.) || (phimax > -Geom::pi()/2.) ){
66  edm::LogError("TkDetLayers") << " something strange going on, please check " << phimin << " " << phimax << " " << phiWin ;
67  }
68  //edm::LogInfo(TkDetLayers) << " Wedge at pi: phi " << phimin << " " << phimax << " " << phiWin
69  // << " " << 2.*Geom::pi()+phiWin << " " ;
70  phiWin += 2.*Geom::pi();
71  phiPos += Geom::pi();
72  }
73 
74  LocalVector localPos( rmed*cos(phiPos), rmed*sin(phiPos), zPos);
75 
76  LogDebug("TkDetLayers") << "localPos in computeBounds: " << localPos << "\n"
77  << "rmin: " << rmin << "\n"
78  << "rmax: " << rmax << "\n"
79  << "zmin: " << zmin << "\n"
80  << "zmax: " << zmax << "\n"
81  << "phiWin: " << phiWin ;
82 
83  return make_pair(new DiskSectorBounds(rmin,rmax,zmin,zmax,phiWin),
84  plane.toGlobal(localPos) );
85 
86  }
87 
88 
90  computeRotation( const vector<const GeomDet*>& dets,
91  const Surface::PositionType& meanPos)
92  {
93  const Plane& plane = dets.front()->surface();
94 
98 
99  GlobalVector planeXAxis = plane.toGlobal( LocalVector( 1, 0, 0));
100  const GlobalPoint& planePosition = plane.position();
101 
102  if(planePosition.x()*planeXAxis.x()+planePosition.y()*planeXAxis.y() > 0.){
103  yAxis = planeXAxis;
104  }else{
105  yAxis = -planeXAxis;
106  }
107 
108  GlobalVector planeZAxis = plane.toGlobal( LocalVector( 0, 0, 1));
109  if(planeZAxis.z()*planePosition.z() > 0.){
110  zAxis = planeZAxis;
111  }else{
112  zAxis = -planeZAxis;
113  }
114 
115  xAxis = yAxis.cross( zAxis);
116 
117  return Surface::RotationType( xAxis, yAxis);
118  }
119 
120 }
121 
123 BladeShapeBuilderFromDet::build( const vector<const GeomDet*>& dets)
124 {
125  // find mean position
127  Vector posSum(0,0,0);
128  for (vector<const GeomDet*>::const_iterator i=dets.begin(); i!=dets.end(); i++) {
129  posSum += (**i).surface().position().basicVector();
130  }
131  Surface::PositionType meanPos( 0.,0.,posSum.z()/float(dets.size()) );
132 
133  // temporary plane - for the computation of bounds
134  Surface::RotationType rotation = computeRotation( dets, meanPos);
135  Plane tmpPlane( meanPos, rotation);
136 
137 
138  auto bo =
139  computeBounds( dets,tmpPlane );
140  GlobalPoint pos = meanPos+bo.second;
141  //edm::LogInfo(TkDetLayers) << "global pos in operator: " << pos ;
142  return new BoundDiskSector( pos, rotation, bo.first);
143 }
#define LogDebug(id)
GlobalPoint toGlobal(const Point2DBase< Scalar, LocalTag > lp) const
Definition: Surface.h:106
Local3DVector LocalVector
Definition: LocalVector.h:12
T perp() const
Definition: PV3DBase.h:72
ROOT::Math::Plane3D::Vector Vector
Definition: EcalHitMaker.cc:29
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
T y() const
Definition: PV3DBase.h:63
Definition: Plane.h:17
float float float z
LocalPoint toLocal(const GlobalPoint &gp) const
Vector3DBase< typename PreciseFloatType< T, U >::Type, FrameTag > cross(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:119
T z() const
Definition: PV3DBase.h:64
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
T min(T a, T b)
Definition: MathUtil.h:58
static BoundDiskSector * build(const std::vector< const GeomDet * > &dets) __attribute__((cold))
static std::vector< GlobalPoint > corners(const Plane &)
Definition: BoundingBox.cc:24
TkRotation< float > RotationType
constexpr double pi()
Definition: Pi.h:31
static uInt32 F(BLOWFISH_CTX *ctx, uInt32 x)
Definition: blowfish.cc:281
T x() const
Definition: PV3DBase.h:62
const PositionType & position() const