test_WireSelector.cxx
Go to the documentation of this file.
1 // test_WireSelector.cxx
2 
3 #ifndef test_WireSelector_CXX
4 #define test_WireSelector_CXX
5 
6 // David Adams
7 // November 2016
8 //
9 // Test the DUNE geometry including channel mapping.
10 //
11 // This test demonstrates how to configure and use the LArSoft Geometry
12 // service outside the art framework. DUNE geometry and geometry helper
13 // service are used.
14 //
15 // Note the geometry service requires the experiment-specific geometry
16 // helper with the channel map also be loaded.
17 //
18 // The functions that set the expected values must be defined externally.
19 
20 #undef NDEBUG
21 
23 #include <string>
24 #include <iostream>
25 #include <fstream>
26 #include <sstream>
27 #include <iomanip>
30 #include "../WireSelector.h"
31 
32 #include "TGraph.h"
33 #include "TH2F.h"
36 
37 using std::string;
38 using std::cout;
39 using std::cerr;
40 using std::endl;
41 using std::ofstream;
42 using std::istringstream;
43 using std::ostringstream;
44 using std::setw;
45 using std::vector;
46 using std::abs;
47 using geo::View_t;
48 using geo::SigType_t;
49 using geo::CryostatID;
50 using geo::TPCID;
51 using geo::PlaneID;
52 using geo::WireID;
54 using readout::ROPID;
55 using geo::Geometry;
56 using geo::CryostatGeo;
57 using geo::TPCGeo;
58 using geo::WireGeo;
60 
61 string toString(const TVector3& xyz, int w =9) {
62  ostringstream sout;
63  sout.precision(3);
64  sout << "("
65  << setw(w) << std::fixed << xyz.x()
66  << ", "
67  << setw(w) << std::fixed << xyz.y()
68  << ", "
69  << setw(w) << std::fixed << xyz.z()
70  << ")";
71  return sout.str();
72 }
73 
74 string toString(const WireSelector::WireInfo& dat, int w =9) {
75  ostringstream sout;
76  sout.precision(3);
77  sout << "("
78  << setw(w) << std::fixed << dat.x
79  << ", "
80  << setw(w) << std::fixed << dat.y
81  << ", "
82  << setw(w) << std::fixed << dat.z
83  << ") "
84  << " ["
85  << setw(5) << dat.channel
86  << "] "
87  << " ("
88  << setw(w) << std::fixed << dat.driftMax
89  << ", "
90  << setw(w) << std::fixed << dat.length
91  << ", "
92  << setw(w) << std::fixed << dat.pitch
93  << ")";
94  return sout.str();
95 }
96 int test_WireSelector(string gname, double wireAngle, double minDrift, unsigned int nShow, int sigopt) {
97  const string myname = "test_WireSelector: ";
98  cout << myname << "Starting test" << endl;
99 #ifdef NDEBUG
100  cout << myname << "NDEBUG must be off." << endl;
101  abort();
102 #endif
103  string line = "-----------------------------";
104 
105  cout << myname << line << endl;
106  cout << myname << " Geometry: " << gname << endl;
107 
109  if ( wireAngle == 100 ) view = geo::kU;
110  if ( wireAngle == 101 ) view = geo::kV;
111  if ( wireAngle == 102 ) view = geo::kZ;
112  bool useView = wireAngle >= 100;
113 
114  cout << myname << line << endl;
115  cout << myname << "Create configuration." << endl;
116  bool useExistingFcl = false;
117  if (useExistingFcl) {
118  ArtServiceHelper::load_services("test_WireSelector.fcl",
120  }
121  else {
122  std::stringstream config;
123  config << "#include \"geometry_dune.fcl\"" << endl;
124  config << "services.Geometry: @local::" + gname << endl;
125  config << "services.ExptGeoHelperInterface: @local::dune_geometry_helper" << endl;
127  }
128 
130 
131  cout << myname << line << endl;
132  cout << myname << "Input arguments:" << endl;
133  cout << myname << " gname: " << gname << endl;
134  cout << myname << " wireAngle: " << wireAngle << endl;
135  cout << myname << " minDrift: " << minDrift << endl;
136  cout << myname << " nShow: " << nShow << endl;
137  cout << myname << " sigopt: " << sigopt << endl;
138 
139  cout << myname << line << endl;
140  cout << myname << "Construct selector." << endl;
142  if ( useView ) {
143  if ( view != geo::kUnknown ) {
144  cout << myname << " Selecting view " << view << endl;
145  ws.selectView(view);
146  }
147  } else {
148  cout << myname << " Selecting wire angle " << wireAngle << endl;
149  ws.selectWireAngle(wireAngle);
150  }
151  if ( minDrift > 0.0 ) {
152  cout << myname << " Selecting min drift " << minDrift << " cm" << endl;
153  ws.selectDrift(minDrift);
154  }
155  assert( ws.geometry() != nullptr );
156 
157  cout << myname << line << endl;
158  cout << myname << "Selector properties:" << endl;
159  cout << myname << " Geometry name: " << pgeo->DetectorName() << endl;
160  cout << myname << " View: " << ws.view() << endl;
161  cout << myname << " Wire angle: " << ws.wireAngle() << endl;
162  cout << myname << " Wire angle tol: " << ws.wireAngleTolerance() << endl;
163  cout << myname << " Drift min: " << ws.driftMin() << endl;
164  cout << myname << " # planes: " << ws.planeIDs().size() << endl;
165 
166  cout << myname << line << endl;
167  cout << myname << "Geometry properties:" << endl;
168  Index ncha = pgeo->Nchannels();
169  cout << myname << " # channels: " << ncha << endl;
170  const double piOver2 = 0.5*acos(-1.0);
171  Index nwirSel = 0;
172  for ( WireSelector::PlaneID pid : ws.planeIDs() ) {
173  const geo::PlaneGeo& gpla = ws.geometry()->Plane(pid);
174  double thtx = gpla.ThetaZ() - piOver2;
175  Index nwir = gpla.Nwires();
176  const geo::TPCGeo& gtpc = pgeo->TPC(pid);
177  double driftSize = gtpc.ActiveWidth();
178  cout << myname << pid
179  << " view=" << gpla.View()
180  << " thtx=" << std::fixed << std::setprecision(3) << thtx
181  << " nwir=" << nwir
182  << " drift distance =" << driftSize << " cm"
183  << endl;
184  TVector3 wdir = gpla.GetWireDirection();
185  TVector3 ndir = gpla.GetNormalDirection();
186  cout << myname << " Wire direction: " << toString(wdir) << endl;
187  cout << myname << " Normal direction: " << toString(ndir) << endl;
188  Index icha1 = ncha;
189  Index icha2 = ncha;
190  Index icha = ncha;
191  string wireText;
192  bool usePlane = useView ? gpla.View() == view : fabs(thtx - wireAngle)<0.001;
193  if ( ! usePlane ) continue;
194  nwirSel += nwir;
195  for ( Index iwir=0; iwir<=nwir; ++iwir ) {
196  bool showWire = iwir < nShow;
197  bool endPlane = iwir == nwir;
198  bool endBlock = endPlane;
199  ostringstream wout;
200  if ( ! endBlock ) {
201  WireID wid(pid, iwir);
202  const WireGeo* pgwir = pgeo->WirePtr(wid);
203  TVector3 xyzWire = pgwir->GetCenter<TVector3>();
204  TVector3 xyz1 = pgwir->GetStart<TVector3>();
205  TVector3 xyz2 = pgwir->GetEnd<TVector3>();
206  icha = pgeo->PlaneWireToChannel(wid);
207  wout << myname << setw(12) << iwir << " [" << setw(4) << icha << "] " << toString(xyzWire);
208  xyzWire.RotateX(thtx);
209  wout << " --> " << toString(xyzWire) << "\n";
210  wout << myname << setw(20) << "" << toString(xyz1);
211  xyz1.RotateX(thtx);
212  wout << " --> " << toString(xyz1) << "\n";
213  wout << myname << setw(20) << "" << toString(xyz2);
214  xyz2.RotateX(thtx);
215  wout << " --> " << toString(xyz2) << "\n";
216  assert( icha < ncha );
217  if ( icha1 == ncha ) {
218  icha1 = icha;
219  icha2 = icha;
220  wireText = showWire ? wout.str() : "";
221  } else if ( icha == icha2 + 1 ) {
222  icha2 = icha;
223  if ( showWire ) wireText += wout.str();
224  } else {
225  endBlock = true;
226  }
227  }
228  if ( endBlock ) {
229  cout << myname << " [" << icha1 << ", " << icha2 << "]" << endl;
230  if ( endPlane ) {
231  icha1 = ncha;
232  icha2 = ncha;
233  } else {
234  icha1 = icha;
235  icha2 = icha;
236  }
237  cout << wireText;
238  wireText = showWire ? wout.str() : "";
239  }
240  }
241  assert( icha1 == ncha );
242  assert( icha2 == ncha );
243  }
244 
245 
246  cout << myname << line << endl;
247  cout << myname << "Check geometry: " << gname << endl;
248  assert( ws.geometry() == pgeo.get() );
249 
250  cout << myname << line << endl;
251  cout << myname << "Check data" << endl;
252  ws.fillData();
253  ws.fillDataMap();
254  cout << myname << " Expected wire count: " << nwirSel << endl;
255  cout << myname << " Reported wire count: " << ws.data().size() << endl;
256  cout << myname << " Mapped wire count: " << ws.dataMap().size() << endl;
257  const WireSelector::WireSummary& wsum = ws.fillWireSummary();
258  assert( ws.data().size() == nwirSel );
259  assert( ws.dataMap().size() == nwirSel );
260  assert( wsum.size() == nwirSel );
261 
262  // Build discriminated adcdata as a vector of x-values for each channel.
263  cout << myname << line << endl;
264  cout << myname << "Build ADC data." << endl;
265  vector<vector<float>> adcdata(ncha);
266  if ( sigopt == 2 ) {
267  cout << myname << " Strumming selected wires." << endl;
268  Index nwir = ws.data().size();
269  float dx = (wsum.xmax-wsum.xmin)/nwir;
270  float xval = wsum.xmin + 0.5*dx;
271  for ( Index iwir=0; iwir<nwir; ++iwir ) {
272  const WireSelector::WireInfo& dat = ws.data()[iwir];
273  adcdata[dat.channel].push_back(xval);
274  xval += dx;
275  }
276  } else if ( sigopt == 1 ) {
277  cout << myname << " Strumming channels." << endl;
278  float dx = (wsum.xmax-wsum.xmin)/ncha;
279  float xval = wsum.xmin + 0.5*dx;
280  for ( Index icha=0; icha<ncha; ++icha ) {
281  adcdata[icha].push_back(xval);
282  xval += dx;
283  }
284  }
285 
286  // Build a graph of z vs. x using adcdata.
287  cout << myname << line << endl;
288  cout << myname << "Build signal data." << endl;
289  vector<float> xsigs;
290  vector<float> zsigs;
291  for ( Index icha=0; icha<ncha; ++icha ) {
292  // Loop over the wires read out by this channel.
293  auto range = ws.dataMap().equal_range(icha);
294  for ( auto ient=range.first; ient!=range.second; ++ient ) {
295  const WireSelector::WireInfo* pdat = ient->second;
296  float zsig = pdat->z;
297  // Loop over the ticks hit for this channel.
298  for ( float xsig : adcdata[icha] ) {
299  // Create a space point for wire and tick.
300  xsigs.push_back(xsig);
301  zsigs.push_back(zsig);
302  }
303  }
304  }
305  Index nsig = xsigs.size();
306  TGraph gsigxz(nsig, &zsigs[0], &xsigs[0]);
307  TGraph gsigzx(nsig, &xsigs[0], &zsigs[0]);
308  cout << myname << " # signals: " << nsig << endl;
309  Index nShowSig = nShow < nsig ? nShow : nsig;
310  for ( Index isig=0; isig<nShowSig; ++isig ) {
311  ostringstream sout;
312  sout.precision(3);
313  sout << myname << setw(10) << std::fixed << xsigs[isig] << "," << setw(9) << std::fixed << zsigs[isig];
314  cout << sout.str() << endl;
315  }
316 
317  bool drawWires = true;
318  LineColors lc;
319  double ticklen = 0.015;
320  Index wpadx = 1200;
321  Index wpady = 1000;
322  if ( drawWires ) {
323  double b = 20;
324  double xpmin = wsum.xmin;
325  double xpmax = wsum.xmin;
326  if ( gname == "protodune_geo" ) {
327  xpmin = -400.0;
328  xpmax = 400.0;
329  }
330  ostringstream ssttl;
331  ssttl << gname << " ";
332  if ( useView ) {
333  ssttl << "view=" << view;
334  } else {
335  ssttl << "wireAngle=" << std::setprecision(3) << std::fixed << wireAngle;
336  }
337  if ( minDrift > 0.0 ) {
338  ssttl << ", TPCs with drift > " << std::setprecision(0) << std::fixed << minDrift << " cm";
339  }
340  string sttlxz = ssttl.str() + "; z [cm]; x [cm]";
341  string sttlzx = ssttl.str() + "; x [cm]; z [cm]";
342  TH2* phaxz = new TH2F("phaxz", sttlxz.c_str(), 1, wsum.zmin-b, wsum.zmax+b, 1, xpmin, xpmax);
343  TH2* phazx = new TH2F("phazx", sttlzx.c_str(), 1, xpmin, xpmax, 1, wsum.zmin-b, wsum.zmax+b);
344  phaxz->SetStats(0);
345  phazx->SetStats(0);
346  TGraph gxz(wsum.size(), &wsum.zWire[0], &wsum.xWire[0]);
347  TGraph gzx(wsum.size(), &wsum.xWire[0], &wsum.zWire[0]);
348  gxz.SetMarkerColor(lc.red());
349  gzx.SetMarkerColor(lc.red());
350  TGraph gxzd(wsum.size(), &wsum.zWire[0], &wsum.xCathode[0]);
351  TGraph gzxd(wsum.size(), &wsum.xCathode[0], &wsum.zWire[0]);
352  gxzd.SetMarkerColor(lc.green());
353  gzxd.SetMarkerColor(lc.green());
354  TPadManipulator padxz(wpadx, wpady);
355  padxz.add(phaxz);
356  padxz.add(&gxz, "P");
357  padxz.add(&gxzd, "P");
358  if ( sigopt ) padxz.add(&gsigxz, "P");
359  padxz.addAxis();
360  padxz.setTickLength(ticklen);
361  padxz.print("test_WireSelector_xz.png");
362  TPadManipulator padzx(wpadx, wpady);
363  padzx.add(phazx);
364  padzx.add(&gzx, "P");
365  padzx.add(&gzxd, "P");
366  if ( sigopt ) padzx.add(&gsigzx, "P");
367  padzx.addAxis();
368  padzx.setTickLength(ticklen);
369  padzx.print("test_WireSelector_zx.png");
370  }
371 
372  cout << myname << line << endl;
373  cout << myname << "Done." << endl;
374  return 0;
375 }
376 
377 //**********************************************************************
378 
379 int main(int argc, const char* argv[]) {
380  string gname = "protodune_geo";
381  double thtx = 102;
382  double minDrift = 0.0;
383  unsigned int nShow = 0;
384  int sigopt = 0;
385  if ( argc > 1 ) {
386  string sarg = argv[1];
387  if ( sarg == "-h" ) {
388  cout << "Usage: " << argv[0] << " [gname] [thtx] [minDrift] [nShow] [sigopt]" << endl;
389  cout << " gname: Geometry name, e.g. protodune_geo" << endl;
390  cout << " thtx: Wire angle, e.g. 0 for vertical" << endl;
391  cout << " >= 100 means use view thtx-100, e.g. 2 for kZ" << endl;
392  cout << " minDrift: select TPC with drift > midDrift cm" << endl;
393  cout << " nShow: # wires, signals to print" << endl;
394  cout << " sigopt: Option to add signals to the displays" << endl;
395  cout << " 0 for no signal" << endl;
396  cout << " 1 for strum of detector channels" << endl;
397  cout << " 2 for strum of selected wires" << endl;
398  return 0;
399  }
400  gname = sarg;
401  if ( argc > 2 ) {
402  istringstream ssarg(argv[2]);
403  ssarg >> thtx;
404  if ( argc > 3 ) {
405  istringstream ssarg(argv[3]);
406  ssarg >> minDrift;
407  if ( argc > 4 ) {
408  istringstream ssarg(argv[4]);
409  ssarg >> nShow;
410  if ( argc > 5 ) {
411  istringstream ssarg(argv[5]);
412  ssarg >> sigopt;
413  }
414  }
415  }
416  }
417  }
418  test_WireSelector(gname, thtx, minDrift, nShow, sigopt);
419  cout << "Tests concluded." << endl;
420  return 0;
421 }
422 
423 //**********************************************************************
424 
425 #endif
void GetStart(double *xyz) const
Definition: WireGeo.h:157
double wireAngleTolerance() const
Definition: WireSelector.h:116
Geometry description of a TPC wireThe wire is a single straight segment on a wire plane...
Definition: WireGeo.h:65
const WireInfoMap & dataMap() const
Definition: WireSelector.h:128
IDparameter< geo::CryostatID > CryostatID
Member type of validated geo::CryostatID parameter.
PlaneGeo const & Plane(unsigned int const p, unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified wire.
double wireAngle() const
Definition: WireSelector.h:115
T * get() const
Definition: ServiceHandle.h:63
int add(unsigned int ipad, TObject *pobj, std::string sopt="", bool replace=false)
enum geo::_plane_proj View_t
Enumerate the possible plane projections.
std::string string
Definition: nybbler.cc:12
Planes which measure V.
Definition: geo_types.h:130
Unknown view.
Definition: geo_types.h:136
std::vector< float > xWire
Definition: WireSelector.h:80
const PlaneIDVector & planeIDs() const
Definition: WireSelector.h:123
static ColorType red()
Definition: LineColors.h:27
static constexpr FileOnPath_t FileOnPath
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
Geometry information for a single TPC.
Definition: TPCGeo.h:38
const WireInfoVector & data() const
Definition: WireSelector.h:127
struct vector vector
Class identifying a set of TPC sharing readout channels.
Definition: readout_types.h:70
Planes which measure Z direction.
Definition: geo_types.h:132
unsigned int Index
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
Vector GetNormalDirection() const
Returns the direction normal to the plane.
Definition: PlaneGeo.h:442
const GeometryCore * geometry() const
Definition: WireSelector.h:106
static void load_services(std::string const &config)
Q_EXPORT QTSManip setprecision(int p)
Definition: qtextstream.h:343
art framework interface to geometry description
static ColorType green()
Definition: LineColors.h:28
unsigned int Nchannels() const
Returns the number of TPC readout channels in the detector.
double ThetaZ() const
Angle of the wires from positive z axis; .
Definition: PlaneGeo.cxx:726
void selectView(View view)
T abs(T value)
Planes which measure U.
Definition: geo_types.h:129
void selectWireAngle(double wireAngle, double tol=0.001)
View_t View() const
Which coordinate does this plane measure.
Definition: PlaneGeo.h:184
int precision() const
Definition: qtextstream.h:259
readout::ROPID ROPID
IDparameter< geo::WireID > WireID
Member type of validated geo::WireID parameter.
std::string DetectorName() const
Returns a string with the name of the detector, as configured.
static Config * config
Definition: config.cpp:1054
enum geo::_plane_sigtype SigType_t
int test_WireSelector(string gname, double wireAngle, double minDrift, unsigned int nShow, int sigopt)
int setTickLength(double len)
Geometry information for a single wire plane.The plane is represented in the geometry by a solid whic...
Definition: PlaneGeo.h:82
void selectDrift(double dmin, double dmax=1.e20)
size_t size
Definition: lodepng.cpp:55
int addAxis(bool flag=true)
std::vector< float > zWire
Definition: WireSelector.h:82
The geometry of one entire detector, as served by art.
Definition: Geometry.h:196
string toString(const TVector3 &xyz, int w=9)
double driftMin() const
Definition: WireSelector.h:119
Q_EXPORT QTSManip setw(int w)
Definition: qtextstream.h:331
geo::View_t View
Definition: WireSelector.h:36
const WireSummary & fillWireSummary()
raw::ChannelID_t PlaneWireToChannel(WireID const &wireid) const
Returns the ID of the TPC channel connected to the specified wire.
int main(int argc, const char *argv[])
const WireInfoVector & fillData()
std::vector< float > xCathode
Definition: WireSelector.h:81
void GetEnd(double *xyz) const
Definition: WireGeo.h:163
double ActiveWidth() const
Width (associated with x coordinate) of active TPC volume [cm].
Definition: TPCGeo.h:97
void line(double t, double *p, double &x, double &y, double &z)
unsigned int Nwires() const
Number of wires in this plane.
Definition: PlaneGeo.h:269
View view() const
Definition: WireSelector.h:112
static bool * b
Definition: config.cpp:1043
detail::Node< FrameID, bool > PlaneID
Definition: CRTID.h:125
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
IDparameter< geo::TPCID > TPCID
Member type of validated geo::TPCID parameter.
void GetCenter(double *xyz, double localz=0.0) const
Fills the world coordinate of a point on the wire.
Definition: WireGeo.cxx:73
readout::TPCsetID APAID
const WireInfoMap & fillDataMap()
unsigned int Index
Definition: WireSelector.h:33
int print(std::string fname, std::string spat="{,}")
QTextStream & endl(QTextStream &s)
Vector GetWireDirection() const
Returns the direction of the wires.
Definition: PlaneGeo.h:513
WireGeo const * WirePtr(geo::WireID const &wireid) const
Returns the specified wire.