CTreeGeometry_module.cc
Go to the documentation of this file.
1 // Dump TPC / Wire Geometry and mapping
2 // Chao Zhang (chao@bnl.gov) 2/7/2018
3 // adapted by Wenqiang Gu (wgu@bnl.gov) 8/30/2020
4 
5 #ifndef CTreeGeometry_module
6 #define CTreeGeometry_module
7 
8 // LArSoft includes
9 //#include "lardata/Utilities/DetectorProperties.h"
11 // #include "Utilities/LArProperties.h"
22 // #include "nusimdata/SimulationBase/MCParticle.h"
23 // #include "nusimdata/SimulationBase/MCTruth.h"
25 #include "lardataobj/RawData/raw.h"
30 
31 
32 // Framework includes
39 // #include "art/Framework/Services/Optional/TFileService.h"
42 #include "canvas/Persistency/Common/FindManyP.h"
45 #include "fhiclcpp/ParameterSet.h"
46 
47 
48 // C++ Includes
49 #include <map>
50 #include <vector>
51 // #include <algorithm>
52 #include <fstream>
53 #include <iostream>
54 #include <iomanip>
55 // #include <string>
56 // #include <sstream>
57 // #include <cmath>
58 
59 // #ifdef __MAKECINT__
60 // #pragma link C++ class vector<vector<int> >+;
61 // #pragma link C++ class vector<vector<float> >+;
62 // #endif
63 
64 using namespace std;
65 
66 namespace {
67 
68 class CTreeGeometry : public art::EDAnalyzer {
69 public:
70 
71  explicit CTreeGeometry(fhicl::ParameterSet const& pset);
72  virtual ~CTreeGeometry();
73 
74  void beginJob();
75  void endJob();
76  void analyze(const art::Event& evt);
77 
78  void reconfigure(fhicl::ParameterSet const& pset);
79  void saveChannelWireMap();
80  void printGeometry();
81 
82 
83 private:
84 
85  // the parameters we'll read from the .fcl
86  bool fSaveChannelWireMap;
87 
89 
90  // Geometry Tree Leafs
91  int fNcryostats;
92  int fNTPC;
93  vector<float> fTPC_x; // TPC length in x
94  vector<float> fTPC_y; // TPC length in y
95  vector<float> fTPC_z; // TPC length in z
96  // int fNplanes; // unused
97  vector<int> fPlane_type; // plane type: 0 == induction, 1 == collection
98  vector<int> fPlane_view; // wire orientation: 0 == U, 1 == V, 2 == X
99  vector<double> fPlane_wirepitch; // wire pitch of each plane
100  vector<double> fPlane_wireangle; // wire angle (to vertical) of each plane
101  vector<int> fPlane_wires; // number of wires in each plane
102  int fNchannels;
103  vector<int> channel_starts; // vector os channel starts on each plane
104  vector<int> channel_ends; // vector os channel starts on each plane
105  //int fNOpDets; // unused
106 
107  // Event Tree Leafs
108  //int fEvent; // unused
109  //int fRun; // unused
110  //int fSubRun; // unused
111 
112  }; // class CTreeGeometry
113 
114 
115 //-----------------------------------------------------------------------
116 CTreeGeometry::CTreeGeometry(fhicl::ParameterSet const& parameterSet)
117  : EDAnalyzer(parameterSet)
118 {
119  reconfigure(parameterSet);
120 }
121 
122 
123 //-----------------------------------------------------------------------
124 CTreeGeometry::~CTreeGeometry()
125 {
126 }
127 
128 
129 //-----------------------------------------------------------------------
131  fSaveChannelWireMap = p.get< bool >("saveChannelWireMap");
132 }
133 
134 
135 //-----------------------------------------------------------------------
137 {
138 
139  fNcryostats = fGeom->Ncryostats(); // 1
140 
141  fNTPC = fGeom->NTPC();
142  for (int i=0; i<fNTPC; i++) {
143  fTPC_x.push_back(fGeom->DetHalfWidth(i)*2);
144  fTPC_y.push_back(fGeom->DetHalfHeight(i)*2);
145  fTPC_z.push_back(fGeom->DetLength(i));
146  }
147 
148  // fNplanes = fGeom->Nplanes();
149  // cout << "#planes: " << fNplanes << endl;
150  // for (int i=0; i<fNplanes; i++) {
151  // fPlane_type.push_back(fGeom->SignalType(geo::PlaneID(0, 0, i)));
152  // fPlane_view.push_back(fGeom->Plane(i).View());
153  // // fPlane_wirepitch[i] = fGeom->WirePitch(fPlane_view[i]); // this doesn't seem to return the correct value!
154  // fPlane_wirepitch.push_back(fGeom->WirePitch(fPlane_view[i], 1, 0)); // this doesn't seem to return the correct value);
155  // fPlane_wireangle.push_back(fGeom->WireAngleToVertical(fGeom->Plane(i).View()));
156  // fPlane_wires.push_back(fGeom->Nwires(i));
157  // }
158 
159  fNchannels = fGeom->Nchannels();
160  // cout << "#channel: " << fNchannels << endl;
161 
162  // Save Channel Map to text file.
163  if (fSaveChannelWireMap) {
164  saveChannelWireMap();
165  }
166 
167  printGeometry();
168 
169 
170 }
171 
172 
173 //-----------------------------------------------------------------------
174 void CTreeGeometry::saveChannelWireMap()
175 {
176  ofstream out;
177  out.open("ChannelWireGeometry.txt");
178  double xyzStart[3];
179  double xyzEnd[3];
180  out << "# " << fGeom->GDMLFile() << "\n";
181  out << "# channel\ttpc\tplane\twire\tsx\tsy\tsz\tex\tey\tez\n";
182  int current_plane = 0;
183  channel_starts.push_back(0);
184  channel_ends.push_back(0);
185  for (int i=0; i<fNchannels; i++) {
186  std::vector<geo::WireID> wireids = fGeom->ChannelToWire(i);
187  int nWires = wireids.size();
188  for (int j=0; j<nWires; j++) {
189  geo::WireID wid = wireids.at(j);
190  int cstat = wid.Cryostat;
191  int tpc = wid.TPC;
192  int plane = wid.Plane;
193  int wire = wid.Wire;
194 
195  int plane_id = plane + tpc*10 + cstat*100;
196  if (plane_id != current_plane) {
197  current_plane = plane_id;
198  channel_starts.push_back(i);
199  channel_ends.push_back(i);
200  }
201  else {
202  channel_ends[channel_ends.size()-1] = i;
203  }
204 
205  fGeom->WireEndPoints(cstat, tpc, plane, wire, xyzStart, xyzEnd);
206 
207  out << i << "\t" << cstat*2+tpc << "\t" << plane << "\t" << wire << "\t";
208  for (int i=0; i<3; i++) {
209  out << xyzStart[i] << "\t";
210  }
211  for (int i=0; i<3; i++) {
212  out << xyzEnd[i] << "\t";
213  }
214  out << "\n";
215  }
216  }
217  out.close();
218 
219 }
220 
221 
222 //-----------------------------------------------------------------------
223 void CTreeGeometry::endJob()
224 {
225 }
226 
227 
228 //-----------------------------------------------------------------------
229 void CTreeGeometry::printGeometry()
230 {
231  cout << "Detector Name: " << fGeom->DetectorName() << endl;
232  cout << "GDML file: " << fGeom->GDMLFile() << endl;
233  // cout << "fNTPC: " << fNTPC << endl;
234  // for (int i=0; i<fNTPC; i++) {
235  // cout << "\tTPC " << i << ": " << fTPC_x[i] << ", " << fTPC_y[i] << ", " << fTPC_z[i] << endl;
236  // }
237  cout << "TPC (Active) Locations: " << endl;
238  for (geo::TPCGeo const& TPC: fGeom->IterateTPCs()) {
239  // get center in world coordinates
240  double origin[3] = {0.};
241  double center[3] = {0.};
242  TPC.LocalToWorld(origin, center);
243  double tpcDim[3] = {TPC.ActiveHalfWidth(), TPC.ActiveHalfHeight(), 0.5*TPC.ActiveLength() };
244  double xmin = center[0] - tpcDim[0];
245  double xmax = center[0] + tpcDim[0];
246  double ymin = center[1] - tpcDim[1];
247  double ymax = center[1] + tpcDim[1];
248  double zmin = center[2] - tpcDim[2];
249  double zmax = center[2] + tpcDim[2];
250  cout << "\t[" << xmin << ", " << xmax << ", " << ymin << ", " << ymax
251  << ", " << zmin << ", " << zmax << "]" << endl;
252  } // for all TPC
253 
254  int size = channel_starts.size();
255  cout << size << " planes: first channels: ";
256  for (int i=0; i<size ; i++) {
257  cout << channel_starts[i] << ", ";
258  }
259  cout << endl;
260  size = channel_ends.size();
261  cout << size << " planes: last channels: ";
262  for (int i=0; i<size ; i++) {
263  cout << channel_ends[i] << ", ";
264  }
265  cout << endl;
266  // for (geo::TPCGeo const& TPC: fGeom->IterateTPCs()) {
267  // for (geo::PlaneGeo const& plane: TPC.IteratePlanes()) {
268  // cout << plane.ID() <<
269  // " channels: " fGeom->PlaneWireToChannel(plane.FirstWire().node())<< endl;
270  // }
271 
272  // }
273 
274  // cout << "fNplanes: " << fNplanes << endl;
275  // for (int i=0; i<fNplanes; i++) {
276  // cout
277  // << "\tplane " << i
278  // << "( type: " << fPlane_type[i]
279  // << ", view: " << fPlane_view[i]
280  // << ", wirepitch: " << fPlane_wirepitch[i]
281  // << ", wire angle: " << fPlane_wireangle[i]
282  // << ", wires: " << fPlane_wires[i]
283  // << ")" << endl;
284  // }
285  cout << "fNchannels: " << fGeom->Nchannels() << endl;
286  cout << "fNOpDet: " << fGeom->NOpDets() << endl;
287  cout << "fAuxDetectors: " << fGeom->NAuxDets() << endl;
288  cout << endl;
289 }
290 
291 //-----------------------------------------------------------------------
293 {
294 
295  // fEvent = event.id().event();
296  // fRun = event.run();
297  // fSubRun = event.subRun();
298 
299  // printEvent();
300 
301 }
302 
303 DEFINE_ART_MODULE(CTreeGeometry)
304 
305 } // namespace
306 
307 #endif // CTreeGeometry_module
def analyze(root, level, gtrees, gbranches, doprint)
Definition: rootstat.py:69
Encapsulate the construction of a single cyostat.
Geometry information for a single TPC.
Definition: TPCGeo.h:38
CryostatID_t Cryostat
Index of cryostat.
Definition: geo_types.h:212
STL namespace.
WireID_t Wire
Index of the wire within its plane.
Definition: geo_types.h:580
art framework interface to geometry description
decltype(auto) constexpr size(T &&obj)
ADL-aware version of std::size.
Definition: StdUtils.h:92
#define DEFINE_ART_MODULE(klass)
Definition: ModuleMacros.h:67
void beginJob()
Definition: Breakpoints.cc:14
virtual void reconfigure(fhicl::ParameterSet const &pset)
T get(std::string const &key) const
Definition: ParameterSet.h:271
p
Definition: test.py:223
PlaneID_t Plane
Index of the plane within its TPC.
Definition: geo_types.h:493
Definition of data types for geometry description.
Encapsulate the geometry of an optical detector.
Declaration of signal hit object.
Encapsulate the construction of a single detector plane.
def center(depos, point)
Definition: depos.py:117
Provides recob::Track data product.
Declaration of basic channel signal object.
TCEvent evt
Definition: DataStructs.cxx:7
TPCID_t TPC
Index of the TPC within its cryostat.
Definition: geo_types.h:406
constexpr Point origin()
Returns a origin position with a point of the specified type.
Definition: geo_vectors.h:227
QTextStream & endl(QTextStream &s)
Event finding and building.