PlaneGeo.cxx
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////////////////
2 // \file PlaneGeo.cxx
3 // \brief
4 // \version $Id: PlaneGeo.cxx,v 1.9 2012-11-28 03:25:00 gsdavies Exp $
5 // \author $Author: gsdavies $
6 /////////////////////////////////////////////////////////////////////////
7 
10 
11 #include <algorithm>
12 #include <cmath>
13 #include <cstdlib>
14 #include <iostream>
15 
16 #include "TGeoNode.h"
17 #include "TGeoVolume.h"
18 #include "TGeoMatrix.h"
19 
20 #include "cetlib_except/exception.h"
21 
22 namespace geo
23 {
24  //..............................................................
25  static bool sort_hori(const CellGeo* c1, const CellGeo* c2)
26  {
27  double xyz1[3], xyz2[3];
28  c1->GetCenter(xyz1);
29  c2->GetCenter(xyz2);
30  return xyz1[1]<xyz2[1];
31  }
32  static bool sort_vert(const CellGeo* c1, const CellGeo* c2)
33  {
34  double xyz1[3], xyz2[3];
35  c1->GetCenter(xyz1);
36  c2->GetCenter(xyz2);
37  return xyz1[0]<xyz2[0];
38  }
39 
40  //......................................................................
41 
42  PlaneGeo::PlaneGeo(std::vector<const TGeoNode*>& n, unsigned int depth)
43  {
44  // build a matrix to take us from the local to the world coordinates
45  // in one step
46  fGeoMatrix = new TGeoHMatrix(*n[0]->GetMatrix());
47  for (unsigned int i=1; i<=depth; ++i) {
48  fGeoMatrix->Multiply(n[i]->GetMatrix());
49  }
50 
51  // Copy fGeoMatrix to Original
52  fGeoMatrixOriginal = new TGeoHMatrix(*fGeoMatrix);
53 
54  // fGeoMatrix is still unmodified
55  fIsOriginal = true;
56 
57  // Determine read-out orientation for the plane. By convention, the
58  // planes, cells, and fibers are constructed such that the z-axis
59  // goes along the long direction with z increasing as one moves
60  // toward the read out end. So, a good way to find the read out end
61  // is to take a small step in the +z direction which moves you
62  // toward the read out side of the detector.
63  static const double step = 10.0; // size of step in cm
64  static const double tol = 0.1*step; // testing tolerance
65  static const double a0[3] = {0,0,0};
66  static const double a1[3] = {0,0,step};
67  double b0[3], b1[3];
68  fGeoMatrix->LocalToMaster(a0,b0);
69  fGeoMatrix->LocalToMaster(a1,b1);
70  b1[0] -= b0[0]; b1[1] -= b0[1]; b1[2] -= b0[2];
71  if (fabs(b1[0]-0.0)<tol && fabs(b1[1]-step)<tol) {
72  fView = kX;
73  fReadout = kTop;
74  }
75  else if (fabs(b1[0]-step)<tol && fabs(b1[1]-0.)<tol) {
76  fView = kY;
77  fReadout = kWest;
78  }
79  else if (fabs(b1[0]+step)<tol && fabs(b1[1]-0.)<tol) {
80  fView = kY;
81  fReadout = kEast;
82  }
83  else {
84  throw cet::exception("BAD_PLANE_VIEW")
85  << "Could not determine view for plane\n"
86  << __FILE__ << ":" << __LINE__ << "\n";
87  }
88  this->FindCells(n, depth);
89 
90  // Make sure horizontals are sorted top to bottom and verticals are
91  // sorted east to west
92  if (fView == kX) {
93  std::sort(fCell.begin(), fCell.end(), sort_vert);
94  }
95  else if (fView == kY) {
96  std::sort(fCell.begin(), fCell.end(), sort_hori);
97  }
98  else abort();
99  }
100 
101  //......................................................................
103  {
104  if(fGeoMatrix) delete fGeoMatrix;
106 
107  for(size_t c = 0; c < fCell.size(); ++c){
108  if(fCell[c]) { delete fCell[c]; fCell[c] = 0; }
109  }
110  fCell.clear();
111  }
112 
113  //......................................................................
114 
115  void PlaneGeo::FindCells(std::vector<const TGeoNode*>& n,
116  unsigned int depth)
117  {
118  if (strncmp(n[depth]->GetName(),"vCell",5)==0) {
119  this->MakeCell(n,depth);
120  return;
121  }
122 
123  // Explore the next layer down
124  unsigned int deeper = depth+1;
125  if (deeper>=n.size()) {
126  throw cet::exception("BAD_NODE")
127  << "Exceeded maximum depth\n"
128  << __FILE__ << ":" << __LINE__ << "\n";
129  }
130  const TGeoVolume* v = n[depth]->GetVolume();
131  int nd = v->GetNdaughters();
132  for (int i=0; i<nd; ++i) {
133  n[deeper] = v->GetNode(i);
134  this->FindCells(n, deeper);
135  }
136  }
137 
138  //......................................................................
139 
140  void PlaneGeo::MakeCell(std::vector<const TGeoNode*>& n,
141  unsigned int depth)
142  {
143  fCell.push_back(new CellGeo(n, depth));
144  }
145 
146  //......................................................................
147 
148  void PlaneGeo::LocalToWorld(const double* plane, double* world) const
149  {
150  fGeoMatrix->LocalToMaster(plane, world);
151  }
152 
153  //......................................................................
154 
155  void PlaneGeo::LocalToWorldVect(const double* plane, double* world) const
156  {
157  fGeoMatrix->LocalToMasterVect(plane, world);
158  }
159 
160  //......................................................................
161 
162  void PlaneGeo::WorldToLocal(const double* world, double* plane) const
163  {
164  fGeoMatrix->MasterToLocal(world, plane);
165  }
166 
167  //......................................................................
168 
169  /// Convert a vector from world frame to the local plane frame
170  /// \param world : 3-D array. Vector in world coordinates; input.
171  /// \param plane : 3-D array. Vector in plane coordinates; plane.
172  void PlaneGeo::WorldToLocalVect(const double* world, double* plane) const
173  {
174  fGeoMatrix->MasterToLocalVect(world,plane);
175  }
176 
177  //......................................................................
178 
179  void PlaneGeo::TranslatePlane(double dx, double dy, double dz,
180  double phi, double theta, double psi)
181  {
182 
183  /// phi, theta and psi are euler angles
184  TGeoRotation *rot = new TGeoRotation("alignmentRot",phi,theta,psi);
185  /// dx, dy and dz are the translations we want
186  TGeoTranslation *trans = new TGeoTranslation("",dx,dy,dz);
187  TGeoCombiTrans *comb = new TGeoCombiTrans(*trans,*rot);
188 
189  this->TranslatePlane(comb);
190 
191  }
192 
193  //......................................................................
194 
195  void PlaneGeo::TranslatePlane(TGeoCombiTrans* comb)
196  {
197 
198  /// Restart from original matrix
199  fGeoMatrix = new TGeoHMatrix(*fGeoMatrixOriginal);
200 
201  /// Apply shift
202  fGeoMatrix->MultiplyLeft(comb->MakeClone());
203 
204  /// Now apply shift to all cells
205  for(unsigned int i=0; i<fCell.size(); i++){
206  fCell.at(i)->TranslateCell(comb);
207  }
208 
209  // fGeoMatrix is now modified
210  fIsOriginal = false;
211 
212  }
213 
214  //......................................................................
215 
217  {
218 
219  /// Restore original matrix
220  fGeoMatrix = new TGeoHMatrix(*fGeoMatrixOriginal);
221 
222  // fGeoMatrix is back to original
223  fIsOriginal = true;
224 
225  }
226 
227 } // end namespace geo
228 ////////////////////////////////////////////////////////////////////////
void GetCenter(double *xyz, double localz=0.0) const
Definition: CellGeo.cxx:159
void FindCells(std::vector< const TGeoNode * > &n, unsigned int depth)
Definition: PlaneGeo.cxx:115
Vertical planes which measure X.
Definition: PlaneGeo.h:28
PlaneGeo(std::vector< const TGeoNode * > &n, unsigned int depth)
Construct a representation of a single plane of the detector.
Definition: PlaneGeo.cxx:42
Vertical modules readout on the top.
Definition: PlaneGeo.h:21
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void WorldToLocal(const double *world, double *plane) const
Transform point from world frame to local plane frame.
Definition: PlaneGeo.cxx:162
View_t fView
Does this plane measure X or Y?
Definition: PlaneGeo.h:88
Horizontal modules readout on the west.
Definition: PlaneGeo.h:22
Horizontal planes which measure Y.
Definition: PlaneGeo.h:29
bool fIsOriginal
Is fGeoMatrix the original transform?
Definition: PlaneGeo.h:86
c2
Definition: demo5.py:33
TGeoHMatrix * fGeoMatrixOriginal
Original plane to world transform.
Definition: PlaneGeo.h:85
static bool sort_vert(const CellGeo *c1, const CellGeo *c2)
Definition: PlaneGeo.cxx:32
Readout_t fReadout
Which end is the readout on?
Definition: PlaneGeo.h:87
std::string GetName(int i)
double dy[NP][NC]
double dx[NP][NC]
double dz[NP][NC]
TGeoHMatrix * fGeoMatrix
Plane to world transform.
Definition: PlaneGeo.h:84
TH1F * a1
Definition: f2_nu.C:476
static bool sort_hori(const CellGeo *c1, const CellGeo *c2)
Definition: PlaneGeo.cxx:25
std::void_t< T > n
void TranslatePlane(double dx, double dy=0, double dz=0, double phi=0, double theta=0, double psi=0)
Definition: PlaneGeo.cxx:179
std::vector< CellGeo * > fCell
List of cells in this plane.
Definition: PlaneGeo.h:89
c1
Definition: demo5.py:24
void LocalToWorldVect(const double *plane, double *world) const
Transform direction vector from local to world.
Definition: PlaneGeo.cxx:155
void MakeCell(std::vector< const TGeoNode * > &n, unsigned int depth)
Definition: PlaneGeo.cxx:140
Encapsulate the cell geometry.
Definition: CellGeo.h:25
Helper for AttenCurve.
Definition: Path.h:10
Horizontal modules readout on the east.
Definition: PlaneGeo.h:23
void RestoreOriginal()
Restore geometry to original.
Definition: PlaneGeo.cxx:216
void LocalToWorld(const double *plane, double *world) const
Transform point from local plane frame to world frame.
Definition: PlaneGeo.cxx:148
void WorldToLocalVect(const double *world, double *plane) const
Transform direction vector from world to local.
Definition: PlaneGeo.cxx:172