CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
volumeHandle.cc
Go to the documentation of this file.
1 // #include "Utilities/Configuration/interface/Architecture.h"
2 
3 /*
4  * See header file for a description of this class.
5  *
6  * $Date: 2013/04/15 16:02:44 $
7  * $Revision: 1.15 $
8  * \author N. Amapane - INFN Torino
9  */
10 
17 
22 
23 #include "CLHEP/Units/GlobalSystemOfUnits.h"
24 
26 
27 #include <string>
28 #include <iterator>
29 #include <iomanip>
30 #include <iostream>
31 #include <boost/lexical_cast.hpp>
32 
33 using namespace SurfaceOrientation;
34 using namespace std;
35 
36 
38  delete refPlane;
39 }
40 
42  : name(fv.logicalPart().name().name()),
43  copyno(fv.copyno()),
44  magVolume(0),
45  masterSector(1),
46  theRN(0.),
47  theRMin(0.),
48  theRMax(0.),
49  refPlane(0),
50  solid(fv.logicalPart().solid()),
51  center_(GlobalPoint(fv.translation().x()/cm,
52  fv.translation().y()/cm,
53  fv.translation().z()/cm)),
54  expand(expand2Pi),
55  isIronFlag(false)
56 {
57  // ASSUMPTION: volume names ends with "_NUM" where NUM is the volume number
58  string volName = name;
59  volName.erase(0,volName.rfind('_')+1);
60  volumeno =boost::lexical_cast<unsigned short>(volName);
61 
62  for (int i=0; i<6; ++i) {
63  isAssigned[i] = false;
64  }
65 
66 
68  cout.precision(7);
69  }
70 
71 
72  referencePlane(fv);
73 
74  if (solid.shape() == ddbox) {
75  buildBox(fv);
76  } else if (solid.shape() == ddtrap) {
77  buildTrap(fv);
78  } else if (solid.shape() == ddcons) {
79  buildCons(fv);
80  } else if (solid.shape() == ddtubs) {
81  buildTubs(fv);
82  } else if (solid.shape() == ddpseudotrap) {
83  buildPseudoTrap(fv);
84  } else if (solid.shape() == ddtrunctubs) {
85  buildTruncTubs(fv);
86  } else {
87  cout << "volumeHandle ctor: Unexpected solid: " << (int) solid.shape() << endl;
88  }
89 
90 
91  // NOTE: Table name and master sector are no longer taken from xml!
92 // DDsvalues_type sv(fv.mergedSpecifics());
93 
94 // { // Extract the name of associated field file.
95 // std::vector<std::string> temp;
96 // std::string pname = "table";
97 // DDValue val(pname);
98 // DDsvalues_type sv(fv.mergedSpecifics());
99 // if (DDfetch(&sv,val)) {
100 // temp = val.strings();
101 // if (temp.size() != 1) {
102 // cout << "*** WARNING: volume has > 1 SpecPar " << pname << endl;
103 // }
104 // magFile = temp[0];
105 
106 // string find="[copyNo]";
107 // std::size_t j;
108 // for ( ; (j = magFile.find(find)) != string::npos ; ) {
109 // stringstream conv;
110 // conv << setfill('0') << setw(2) << copyno;
111 // string repl;
112 // conv >> repl;
113 // magFile.replace(j, find.length(), repl);
114 // }
115 
116 // } else {
117 // cout << "*** WARNING: volume does not have a SpecPar " << pname << endl;
118 // cout << " DDsvalues_type: " << fv.mergedSpecifics() << endl;
119 // }
120 // }
121 
122 // { // Extract the number of the master sector.
123 // std::vector<double> temp;
124 // const std::string pname = "masterSector";
125 // DDValue val(pname);
126 // if (DDfetch(&sv,val)) {
127 // temp = val.doubles();
128 // if (temp.size() != 1) {
129 // cout << "*** WARNING: volume has > 1 SpecPar " << pname << endl;
130 // }
131 // masterSector = int(temp[0]+.5);
132 // } else {
133 // if (MagGeoBuilderFromDDD::debug) {
134 // cout << "Volume does not have a SpecPar " << pname
135 // << " using: " << copyno << endl;
136 // cout << " DDsvalues_type: " << fv.mergedSpecifics() << endl;
137 // }
138 // masterSector = copyno;
139 // }
140 // }
141 
142  // Get material for this volume
143  if (fv.logicalPart().material().name().name() == "Iron") isIronFlag=true;
144 
145 
147  cout << " RMin = " << theRMin <<endl;
148  cout << " RMax = " << theRMax <<endl;
149 
150  if (theRMin < 0 || theRN < theRMin || theRMax < theRN)
151  cout << "*** WARNING: wrong RMin/RN/RMax , shape: " << (int) shape() << endl;
152 
153  cout << "Summary: " << name << " " << copyno
154  << " Shape= " << (int) shape()
155  << " trasl " << center()
156  << " R " << center().perp()
157  << " phi " << center().phi()
158  << " magFile " << magFile
159  << " Material= " << fv.logicalPart().material().name()
160  << " isIron= " << isIronFlag
161  << " masterSector= " << masterSector << std::endl;
162 
163  cout << " Orientation of surfaces:";
164  std::string sideName[3] = {"positiveSide", "negativeSide", "onSurface"};
165  for (int i=0; i<6; ++i) {
166  cout << " " << i << ":" << sideName[surfaces[i]->side(center_,0.3)];
167  }
168  cout << endl;
169  }
170 }
171 
172 
174  return center_;
175 }
176 
178  // The refPlane is the "main plane" for the solid. It corresponds to the
179  // x,y plane in the DDD local frame, and defines a frame where the local
180  // coordinates are the same as in DDD.
181  // In the geometry version 85l_030919, this plane is normal to the
182  // beam line for all volumes but pseudotraps, so that global R is along Y,
183  // global phi is along -X and global Z along Z:
184  //
185  // Global(for vol at pi/2) Local
186  // +R (+Y) +Y
187  // +phi(-X) -X
188  // +Z +Z
189  //
190  // For pseudotraps the refPlane is parallel to beam line and global R is
191  // along Z, global phi is along +-X and and global Z along Y:
192  //
193  // Global(for vol at pi/2) Local
194  // +R (+Y) +Z
195  // +phi(-X) +X
196  // +Z +Y
197  //
198  // Note that the frame is centered in the DDD volume center, which is
199  // inside the volume for DDD boxes and (pesudo)trapezoids, on the beam line
200  // for tubs, cons and trunctubs.
201 
202  // In geometry version 1103l, trapezoids have X and Z in the opposite direction
203  // than the above. Boxes are either oriented as described above or in some case
204  // have opposite direction for Y and X.
205 
206  // The global position
207  Surface::PositionType & posResult = center_;
208 
209  // The reference plane rotation
210  DD3Vector x, y, z;
211  fv.rotation().GetComponents(x,y,z);
213  if (x.Cross(y).Dot(z) < 0.5) {
214  cout << "*** WARNING: Rotation is not RH "<< endl;
215  }
216  }
217 
218  // The global rotation
220  rotResult(float(x.X()),float(x.Y()),float(x.Z()),
221  float(y.X()),float(y.Y()),float(y.Z()),
222  float(z.X()),float(z.Y()),float(z.Z()));
223 
224  refPlane = new GloballyPositioned<float>(posResult, rotResult);
225 
226  // Check correct orientation
228 
229  cout << "Refplane pos " << refPlane->position() << endl;
230 
231  // See comments above for the conventions for orientation.
232  LocalVector globalZdir(0.,0.,1.); // Local direction of the axis along global Z
233  if (solid.shape() == ddpseudotrap) {
234  globalZdir = LocalVector(0.,1.,0.);
235  }
236  if (refPlane->toGlobal(globalZdir).z()<0.) {
237  globalZdir=-globalZdir;
238  }
239 
240  float chk = refPlane->toGlobal(globalZdir).dot(GlobalVector(0,0,1));
241  if (chk < .999) cout << "*** WARNING RefPlane check failed!***"
242  << chk << endl;
243  }
244 }
245 
246 
247 
249  double deltaPhi,
250  double zhalf,
251  double rCentr) {
252  // This is 100% equal for cons and tubs!!!
253 
254  GlobalVector planeXAxis = refPlane->toGlobal(LocalVector( 1, 0, 0));
255  GlobalVector planeYAxis = refPlane->toGlobal(LocalVector( 0, 1, 0));
256  GlobalVector planeZAxis = refPlane->toGlobal(LocalVector( 0, 0, 1));
257 
258  // Local Y axis of the faces at +-phi.
259  GlobalVector y_phiplus = refPlane->toGlobal(LocalVector(cos(startPhi+deltaPhi),
260  sin(startPhi+deltaPhi),0.));
261  GlobalVector y_phiminus = refPlane->toGlobal(LocalVector(cos(startPhi),
262  sin(startPhi),0.));
263 
264  Surface::RotationType rot_Z(planeXAxis,planeYAxis);
265  Surface::RotationType rot_phiplus(planeZAxis, y_phiplus);
266  Surface::RotationType rot_phiminus(planeZAxis, y_phiminus);
267 
268  GlobalPoint pos_zplus(center_.x(),center_.y(),center_.z()+zhalf);
269  GlobalPoint pos_zminus(center_.x(),center_.y(),center_.z()-zhalf);
270  // BEWARE: in this case, the origin for phiplus,phiminus surfaces is
271  // at radius R and not on a plane passing by center_ orthogonal to the radius.
272  GlobalPoint pos_phiplus(refPlane->toGlobal(LocalPoint(rCentr*cos(startPhi+deltaPhi),rCentr*sin(startPhi+deltaPhi),0.)));
273  GlobalPoint pos_phiminus(refPlane->toGlobal(LocalPoint(rCentr*cos(startPhi),
274  rCentr*sin(startPhi),
275  0.)));
276  surfaces[zplus] = new Plane(pos_zplus, rot_Z);
277  surfaces[zminus] = new Plane(pos_zminus, rot_Z);
278  surfaces[phiplus] = new Plane(pos_phiplus, rot_phiplus);
279  surfaces[phiminus] = new Plane(pos_phiminus, rot_phiminus);
280 
282  cout << "Actual Center at: " << center_ << " R " << center_.perp()
283  << " phi " << center_.phi() << endl;
284  cout << "RN " << theRN << endl;
285 
286  cout << "pos_zplus " << pos_zplus << " "
287  << pos_zplus.perp() << " " << pos_zplus.phi() << endl
288  << "pos_zminus " << pos_zminus << " "
289  << pos_zminus.perp() << " " << pos_zminus.phi() << endl
290  << "pos_phiplus " << pos_phiplus << " "
291  << pos_phiplus.perp() << " " << pos_phiplus.phi() <<endl
292  << "pos_phiminus " << pos_phiminus << " "
293  << pos_phiminus.perp() << " " << pos_phiminus.phi() <<endl;
294 
295  cout << "y_phiplus " << y_phiplus << endl;
296  cout << "y_phiminus " << y_phiminus << endl;
297 
298  cout << "rot_Z " << surfaces[zplus]->toGlobal(LocalVector(0.,0.,1.)) << endl
299  << "rot_phi+ " << surfaces[phiplus]->toGlobal(LocalVector(0.,0.,1.))
300  << " phi " << surfaces[phiplus]->toGlobal(LocalVector(0.,0.,1.)).phi()
301  << endl
302  << "rot_phi- " << surfaces[phiminus]->toGlobal(LocalVector(0.,0.,1.))
303  << " phi " << surfaces[phiminus]->toGlobal(LocalVector(0.,0.,1.)).phi()
304  << endl;
305  }
306 
307 // // Check ordering.
309  if (pos_zplus.z() < pos_zminus.z()) {
310  cout << "*** WARNING: pos_zplus < pos_zminus " << endl;
311  }
312  if (Geom::Phi<float>(pos_phiplus.phi()-pos_phiminus.phi()) < 0. ) {
313  cout << "*** WARNING: pos_phiplus < pos_phiminus " << endl;
314  }
315  }
316 }
317 
318 
319 
320 bool MagGeoBuilderFromDDD::volumeHandle::sameSurface(const Surface & s1, Sides which_side, float tolerance)
321 {
322  //Check for null comparison
323  if (&s1==(surfaces[which_side]).get()){
324  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: OK (same ptr)" << endl;
325  return true;
326  }
327 
328  const float maxtilt = 0.999;
329 
330  const Surface & s2 = *(surfaces[which_side]);
331  // Try with a plane.
332  const Plane * p1 = dynamic_cast<const Plane*>(&s1);
333  if (p1!=0) {
334  const Plane * p2 = dynamic_cast<const Plane*>(&s2);
335  if (p2==0) {
336  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: different types" << endl;
337  return false;
338  }
339 
340  if ( (fabs(p1->normalVector().dot(p2->normalVector())) > maxtilt)
341  && (fabs((p1->toLocal(p2->position())).z()) < tolerance) ) {
342  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: OK "
343  << fabs(p1->normalVector().dot(p2->normalVector()))
344  << " " << fabs((p1->toLocal(p2->position())).z()) << endl;
345  return true;
346  } else{
347  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: not the same: "
348  << p1->normalVector() << p1->position() << endl
349  << " "
350  << p2->normalVector() << p2->position() << endl
351  << fabs(p1->normalVector().dot(p2->normalVector()))
352  << " " << (p1->toLocal(p2->position())).z()<< endl;
353  return false;
354  }
355  }
356 
357  // Try with a cylinder.
358  const Cylinder * cy1 = dynamic_cast<const Cylinder*>(&s1);
359  if (cy1!=0) {
360  const Cylinder * cy2 = dynamic_cast<const Cylinder*>(&s2);
361  if (cy2==0) {
362  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: different types" << endl;
363  return false;
364  }
365  // Assume axis is the same!
366  if (fabs(cy1->radius() - cy2->radius()) < tolerance) {
367  return true;
368  } else {
369  return false;
370  }
371  }
372 
373  // Try with a cone.
374  const Cone * co1 = dynamic_cast<const Cone*>(&s1);
375  if (co1!=0) {
376  const Cone * co2 = dynamic_cast<const Cone*>(&s2);
377  if (co2==0) {
378  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: different types" << endl;
379  return false;
380  }
381  // FIXME
382  if (fabs(co1->openingAngle()-co2->openingAngle()) < maxtilt
383  && (co1->vertex()-co2->vertex()).mag() < tolerance) {
384  return true;
385  } else {
386  return false;
387  }
388  }
389 
390  if (MagGeoBuilderFromDDD::debug) cout << " sameSurface: unknown surfaces..." << endl;
391  return false;
392 }
393 
394 
395 
397 {
398  //Check for null assignment
399  if (&s1==(surfaces[which_side]).get()){
400  isAssigned[which_side] = true;
401  return true;
402  }
403 
404  if (!sameSurface(s1,which_side)){
405  cout << "***ERROR: setSurface: trying to assign a surface that does not match destination surface. Skipping." << endl;
406  const Surface & s2 = *(surfaces[which_side]);
407  //FIXME: Just planes for the time being!!!
408  const Plane * p1 = dynamic_cast<const Plane*>(&s1);
409  const Plane * p2 = dynamic_cast<const Plane*>(&s2);
410  if (p1!=0 && p2 !=0)
411  cout << p1->normalVector() << p1->position() << endl
412  << p2->normalVector() << p2->position() << endl;
413  return false;
414  }
415 
416 
417  if (isAssigned[which_side]) {
418  if (&s1!=(surfaces[which_side]).get()){
419  cout << "*** WARNING volumeHandle::setSurface: trying to reassign a surface to a different surface instance" << endl;
420  return false;
421  }
422  } else {
423  surfaces[which_side] = &s1;
424  isAssigned[which_side] = true;
425  if (MagGeoBuilderFromDDD::debug) cout << " Volume " << name << " # " << copyno << " Assigned: " << (int) which_side << endl;
426  return true;
427  }
428 
429  return false; // let the compiler be happy
430 }
431 
432 
433 
434 const Surface &
436  return *(surfaces[which_side]);
437 }
438 
439 
440 
441 const Surface &
443  assert(which_side >=0 && which_side <6);
444  return *(surfaces[which_side]);
445 }
446 
447 
448 std::vector<VolumeSide>
450  std::vector<VolumeSide> result;
451  for (int i=0; i<6; ++i){
452  // If this is just a master volume out of wich a 2pi volume
453  // should be built (e.g. central cylinder), skip the phi boundaries.
454  if (expand && (i==phiplus || i==phiminus)) continue;
455 
456  // FIXME: Skip null inner degenerate cylindrical surface
457  if (solid.shape() == ddtubs && i == SurfaceOrientation::inner && theRMin < 0.001) continue;
458 
459  ReferenceCountingPointer<Surface> s = const_cast<Surface*> (surfaces[i].get());
460  result.push_back(VolumeSide(s, GlobalFace(i),
461  surfaces[i]->side(center_,0.3)));
462  }
463  return result;
464 }
465 
466 void MagGeoBuilderFromDDD::volumeHandle::printUniqueNames(handles::const_iterator begin, handles::const_iterator end ) {
467  std::vector<std::string> names;
468  for (handles::const_iterator i = begin;
469  i != end; ++i){
470  names.push_back((*i)->name);
471  }
472 
473  sort(names.begin(),names.end());
474  std::vector<std::string>::iterator i = unique(names.begin(),names.end());
475  int nvols = int(i - names.begin());
476  cout << nvols << " ";
477  copy(names.begin(), i, ostream_iterator<std::string>(cout, " "));
478 
479  cout << endl;
480 }
481 
482 
483 #include "MagneticField/GeomBuilder/src/buildBox.icc"
484 #include "MagneticField/GeomBuilder/src/buildTrap.icc"
485 #include "MagneticField/GeomBuilder/src/buildTubs.icc"
486 #include "MagneticField/GeomBuilder/src/buildCons.icc"
487 #include "MagneticField/GeomBuilder/src/buildPseudoTrap.icc"
488 #include "MagneticField/GeomBuilder/src/buildTruncTubs.icc"
unsigned short copyno
copy number
Definition: volumeHandle.h:64
int i
Definition: DBlmapReader.cc:9
void buildBox(const DDExpandedView &fv)
const DDRotationMatrix & rotation() const
The absolute rotation of the current node.
Local3DVector LocalVector
Definition: LocalVector.h:12
void buildPhiZSurf(double startPhi, double deltaPhi, double zhalf, double rCentr)
const N & name() const
Definition: DDBase.h:82
T perp() const
Definition: PV3DBase.h:72
static const HistoName names[]
T mag() const
The vector magnitude. Equivalent to sqrt(vec.mag2())
void buildTubs(const DDExpandedView &fv)
GlobalVector normalVector() const
Definition: Plane.h:45
std::vector< VolumeSide > sides() const
The surfaces and they orientation, as required to build a MagVolume.
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Geom::Phi< T > phi() const
Definition: PV3DBase.h:69
void buildPseudoTrap(const DDExpandedView &fv)
const GlobalPoint & center() const
Return the center of the volume.
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
static void printUniqueNames(handles::const_iterator begin, handles::const_iterator end)
Just for debugging...
std::string name
Name of the volume.
Definition: volumeHandle.h:58
DDSolidShape shape() const
Shape of the solid.
Definition: volumeHandle.h:99
Definition: Plane.h:17
float float float z
tuple s2
Definition: indexGen.py:106
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< double > > DD3Vector
A DD Translation is currently implemented with Root Vector3D.
Definition: DDTranslation.h:6
LocalPoint toLocal(const GlobalPoint &gp) const
T z() const
Definition: PV3DBase.h:64
tuple result
Definition: query.py:137
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
std::string magFile
Name of magnetic field table file.
Definition: volumeHandle.h:60
DDSolidShape shape(void) const
The type of the solid.
Definition: DDSolid.cc:144
void buildCons(const DDExpandedView &fv)
#define end
Definition: vmac.h:38
void referencePlane(const DDExpandedView &fv)
double p2[4]
Definition: TauolaWrapper.h:90
volumeHandle(const DDExpandedView &fv, bool expand2Pi=false)
Definition: volumeHandle.cc:41
unsigned short volumeno
volume number
Definition: volumeHandle.h:62
const Surface & surface(int which_side) const
Get the current surface on specified side.
#define begin
Definition: vmac.h:31
double p1[4]
Definition: TauolaWrapper.h:89
void buildTrap(const DDExpandedView &fv)
int masterSector
The sector for which an interpolator for this class of volumes should be built.
Definition: volumeHandle.h:113
tuple cout
Definition: gather_cfg.py:121
const DDLogicalPart & logicalPart() const
The logical-part of the current node in the expanded-view.
Definition: DDAxes.h:10
const PositionType & position() const
Provides an exploded view of the detector (tree-view)
bool setSurface(const Surface &s1, Sides which_side)
Assign a shared surface perorming sanity checks.
bool sameSurface(const Surface &s1, Sides which_side, float tolerance=0.01)
Find out if two surfaces are the same physical surface.
const std::string & name() const
Returns the name.
Definition: DDName.cc:87
const DDMaterial & material(void) const
Returns a reference object of the material this LogicalPart is made of.
Global3DVector GlobalVector
Definition: GlobalVector.h:10
void buildTruncTubs(const DDExpandedView &fv)