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