test_GeometryDune10kt.cxx
Go to the documentation of this file.
1 // test_GeometryDune10kt.cxx
2 
3 // David Adams
4 // November 2016
5 //
6 // Test the DUNE 10kt geometry including channel mapping.
7 
8 #include "test_GeometryDune.cxx"
9 
10 //**********************************************************************
11 
12 // Set expected values.
13 
14 void setExpectedValues(ExpectedValues& ev) {
15  ev.gname ="dune10kt_workspace_geo";
16  ev.fullname = "dune10kt_v1_workspace";
17  // Geometry counts.
18  ev.ncry = 1;
19  ev.ntpc = 8;
20  ev.npla = 3;
21  ev.napa = 4;
22  ev.nrop = 4;
23  // Signal type and view for each TPC plane.
24  resize(ev.sigType, ev.ncry, ev.ntpc, ev.npla, geo::kMysteryType);
25  resize(ev.view, 1, 8, 3, geo::kUnknown);
26  for ( Index icry=0; icry<ev.ncry; ++icry ) {
27  for ( Index itpc=0; itpc<ev.ntpc; ++itpc ) {
28  for ( Index ipla=0; ipla<ev.npla; ++ipla ) {
29  ev.sigType[icry][itpc][ipla] = (ipla < 2) ? geo::kInduction : geo::kCollection;
30  }
31  ev.view[icry][itpc][0] = geo::kU;
32  ev.view[icry][itpc][1] = geo::kV;
33  ev.view[icry][itpc][2] = geo::kZ;
34  }
35  }
36  // Wire count for each TPC plane.
37  resize(ev.nwirPerPlane, ev.ntpc, ev.npla, 0);
38  ev.nwirPerPlane[0][0] = 1149;
39  ev.nwirPerPlane[0][1] = 1148;
40  ev.nwirPerPlane[0][2] = 480;
41  for ( Index itpc=1; itpc<ev.ntpc; ++itpc ) {
42  for ( Index ipla=0; ipla<ev.npla; ++ipla ) ev.nwirPerPlane[itpc][ipla] = ev.nwirPerPlane[0][ipla];
43  }
44  // Channel count per ROP.
45  resize(ev.nchaPerRop, ev.nrop, 0);
46  ev.nchaPerRop[0] = 800;
47  ev.nchaPerRop[1] = 800;
48  ev.nchaPerRop[2] = 480;
49  ev.nchaPerRop[3] = 480;
50  for ( Index irop=0; irop<ev.nrop; ++irop ) ev.nchaPerApa += ev.nchaPerRop[irop];
51  ev.nchatot = ev.napa*ev.nchaPerApa;
52  resize(ev.chacry, ev.nchatot, InvalidIndex);
53  resize(ev.chaapa, ev.nchatot, InvalidIndex);
54  resize(ev.charop, ev.nchatot, InvalidIndex);
55  Index icha = 0;
56  for ( Index icry=0; icry<ev.ncry; ++icry ) {
57  for ( Index iapa=0; iapa<ev.napa; ++iapa ) {
58  for ( Index irop=0; irop<ev.nrop; ++irop ) {
59  for ( Index kcha=0; kcha<ev.nchaPerRop[irop]; ++kcha ) {
60  ev.chacry[icha] = icry;
61  ev.chaapa[icha] = iapa;
62  ev.charop[icha] = irop;
63  ++icha;
64  }
65  }
66  }
67  }
68  resize(ev.firstchan, ev.ncry, ev.napa, ev.nrop, raw::InvalidChannelID);
69  Index chan = 0;
70  for ( Index icry=0; icry<ev.ncry; ++icry ) {
71  for ( Index iapa=0; iapa<ev.napa; ++iapa ) {
72  for ( Index irop=0; irop<ev.nrop; ++irop ) {
73  ev.firstchan[icry][iapa][irop] = chan;
74  chan += ev.nchaPerRop[irop];
75  }
76  }
77  }
78  // Space points.
79  #include "setWorkspaceSpacePoints.dat"
80  // Optical detectors.
81  ev.nopdet = 10*ev.napa;
82  ev.nopdetcha.resize(ev.nopdet, 1);
83  ev.nopcha = ev.nopdet;
84  resize(ev.opdetcha, ev.nopdet, 1, 0);
85  for ( Index iopt=0; iopt<ev.nopdet; ++iopt ) {
86  Index icha = iopt;
87  for ( Index ioch=0; ioch<ev.nopdetcha[iopt]; ++ ioch ) {
88  ev.opdetcha[iopt][ioch] = icha++;
89  }
90  }
91 }
92 
93 //**********************************************************************
94 
96  const string myname = "setExpectedValuesSpacePoints: ";
97  string ofname = "setExpectedValuesSpacePoints.dat";
98  cout << myname << "Writing " << ofname << endl;
99  const CryostatGeo& crygeo = pgeo->Cryostat(0);
100  double origin[3] = {0.0};
101  double crypos[3] = {0.0};
102  crygeo.LocalToWorld(origin, crypos);
103  double cxlo = crypos[0] - crygeo.HalfWidth();
104  double cxhi = crypos[0] + crygeo.HalfWidth();
105  double cylo = crypos[1] - crygeo.HalfHeight();
106  double cyhi = crypos[1] + crygeo.HalfHeight();
107  double czlo = crypos[2] - 0.5*crygeo.Length();
108  double czhi = crypos[2] + 0.5*crygeo.Length();
109  cout << "Cryostat limits: "
110  << "(" << cxlo << ", " << cylo << ", " << czlo << "), "
111  << "(" << cxhi << ", " << cyhi << ", " << czhi << ")" << endl;
112  vector<double> zfs = {0.2, 0.3, 0.4, 0.5, 0.6, 0.7 };
113  vector<double> yfs = {0.2, 0.5, 0.8};
114  vector<double> xfs = {0.2, 0.3, 0.6 };
115  ofstream fout(ofname.c_str());
116  for ( double zf : zfs ) {
117  double z = czlo + zf*(czhi-czlo);
118  for ( double yf : yfs ) {
119  double y = cylo + yf*(cyhi-cylo);
120  for ( double xf : xfs ) {
121  double x = cxlo + xf*(cxhi-cxlo);
122  double xyz[3] = {x, y, z};
123  cout << " (" << x << ", " << y << ", " << z << ")" << endl;
124  TPCID tpcid = pgeo->FindTPCAtPosition(xyz);
125  unsigned int itpc = tpcid.TPC;
126  if ( itpc == TPCID::InvalidID ) continue;
127  const TPCGeo& tpcgeo = pgeo->TPC(tpcid);
128  unsigned int npla = tpcgeo.Nplanes();
129  fout << " ev.posXyz.push_back(SpacePoint(" << x << ", " << y << ", " << z << "));" << endl;
130  for ( unsigned int ipla=0; ipla<npla; ++ipla ) {
131  PlaneID plaid(tpcid, ipla);
132  double xwire = pgeo->WireCoordinate(x, y, plaid);
133  fout << " ev.posTpc.push_back(" << itpc << ");" << endl;
134  fout << " ev.posPla.push_back(" << ipla << ");" << endl;
135  fout << " ev.posWco.push_back(" << xwire << ");" << endl;
136  } // end loop over planes
137  } // end loop over x
138  } // end loop over y
139  } // end loop over z
140 }
141 
142 //**********************************************************************
143 
geo::Length_t WireCoordinate(double YPos, double ZPos, geo::PlaneID const &planeid) const
Returns the index of the nearest wire to the specified position.
Who knows?
Definition: geo_types.h:147
Planes which measure V.
Definition: geo_types.h:130
Unknown view.
Definition: geo_types.h:136
unsigned int Nplanes() const
Number of planes in this tpc.
Definition: TPCGeo.h:165
The data type to uniquely identify a Plane.
Definition: geo_types.h:472
Geometry information for a single TPC.
Definition: TPCGeo.h:38
void setExpectedValuesSpacePoints(Geometry *pgeo)
Planes which measure Z direction.
Definition: geo_types.h:132
unsigned int Index
Geometry information for a single cryostat.
Definition: CryostatGeo.h:43
void resize(Vector< T > &vec1, Index n1, const V &val)
const Index InvalidIndex
geo::TPCID FindTPCAtPosition(double const worldLoc[3]) const
Returns the ID of the TPC at specified location.
Planes which measure U.
Definition: geo_types.h:129
constexpr ChannelID_t InvalidChannelID
ID of an invalid channel.
Definition: RawTypes.h:32
Signal from induction planes.
Definition: geo_types.h:145
double HalfWidth() const
Half width of the cryostat [cm].
CryostatGeo const & Cryostat(geo::CryostatID const &cryoid) const
Returns the specified cryostat.
The geometry of one entire detector, as served by art.
Definition: Geometry.h:196
The data type to uniquely identify a TPC.
Definition: geo_types.h:386
double HalfHeight() const
Half height of the cryostat [cm].
void LocalToWorld(const double *cryo, double *world) const
Transform point from local cryostat frame to world frame.
Definition: CryostatGeo.h:387
list x
Definition: train.py:276
TPCGeo const & TPC(unsigned int const tpc=0, unsigned int const cstat=0) const
Returns the specified TPC.
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
double Length() const
Length of the cryostat [cm].
Definition: CryostatGeo.h:107
QTextStream & endl(QTextStream &s)
Signal from collection planes.
Definition: geo_types.h:146
void setExpectedValues(ExpectedValues &ev)