Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Accommodate breaking changes in geometry #28

Merged
merged 1 commit into from
Dec 21, 2022
Merged

Accommodate breaking changes in geometry #28

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
67 changes: 37 additions & 30 deletions dunesim/DetSim/Service/Dune35tSimChannelExtractService_service.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -234,27 +234,29 @@ void Dune35tSimChannelExtractService::init() {
// wire running over the comb are beyond the endpoints. So we need to extrapolate.

double xyzbeg[3],xyzend[3];
int lastwire = 0;
unsigned int lastwire = 0;

// Numbers in comments are from Geometry V3 for debugging purposes.

// APA 0

geo->WireEndPoints(0,0,0,0,xyzbeg,xyzend); // first U wire in TPC 0.
constexpr geo::TPCID apa0{0,0};
constexpr geo::PlaneID apa0_plane0{apa0,0};
geo->WireEndPoints(geo::WireID{apa0_plane0,0},xyzbeg,xyzend); // first U wire in TPC 0.
zcomb2 = xyzbeg[2]; // 0.0
ycomb5 = xyzend[1]; // 113.142

lastwire = geo->Nwires(0,0,0)-1; // 358 in v3
geo->WireEndPoints(0,0,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 0.
lastwire = geo->Nwires(apa0_plane0)-1; // 358 in v3
geo->WireEndPoints(geo::WireID{apa0_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 0.
zcomb5 = xyzend[2]; // 50.8929
ycomb2 = xyzbeg[1]; // -82.9389

geo->WireEndPoints(0,0,1,0,xyzbeg,xyzend); // first V wire in TPC 0.
constexpr geo::PlaneID apa0_plane1{apa0, 1};
geo->WireEndPoints(geo::WireID{apa0_plane1,0},xyzbeg,xyzend); // first V wire in TPC 0.
zcomb4 = xyzend[2]; // 50.5774
ycomb4 = xyzbeg[1]; // 113.142

lastwire = geo->Nwires(1,0,0)-1; // 344 in v3
geo->WireEndPoints(0,0,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 0.
lastwire = geo->Nwires(apa0_plane1)-1; // 344 in v3
geo->WireEndPoints(geo::WireID{apa0_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 0.
zcomb3 = xyzbeg[2]; // 0.3155
ycomb3 = xyzend[1]; // -82.6234

Expand All @@ -270,24 +272,25 @@ void Dune35tSimChannelExtractService::init() {
ycomb1 = ycomb2 - (ycomb3 - ycomb2);
ycomb6 = ycomb5 + (ycomb5 - ycomb4);


// APA 1

geo->WireEndPoints(0,2,0,0,xyzbeg,xyzend); // first U wire in TPC 2.
constexpr geo::TPCID apa1{0, 2};
constexpr geo::PlaneID apa1_plane0{apa1,0};
geo->WireEndPoints(geo::WireID{apa1_plane0,0},xyzbeg,xyzend); // first U wire in TPC 2.
zcomb11 = xyzend[2]; // 102.817
ycomb8 = xyzbeg[1]; // -85.221

lastwire = geo->Nwires(0,2,0)-1; // 194 in v3
geo->WireEndPoints(0,2,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 2.
lastwire = geo->Nwires(apa1_plane0)-1; // 194 in v3
geo->WireEndPoints(geo::WireID{apa1_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 2.
zcomb8 = xyzbeg[2]; // 51.924
ycomb11 = xyzend[1]; // -0.831

geo->WireEndPoints(0,2,1,0,xyzbeg,xyzend); // first V wire in TPC 2.
constexpr geo::PlaneID apa1_plane1{apa1, 1};
geo->WireEndPoints(geo::WireID{apa1_plane1,0},xyzbeg,xyzend); // first V wire in TPC 2.
zcomb9 = xyzbeg[2]; // 52.2395
ycomb9 = xyzend[1]; // -85.222

lastwire = geo->Nwires(1,2,0)-1; // 188 in v3
geo->WireEndPoints(0,2,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 2.
lastwire = geo->Nwires(apa1_plane1)-1; // 188 in v3
geo->WireEndPoints(geo::WireID{apa1_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 2.
zcomb10 = xyzend[2]; // 102.501
ycomb10 = xyzbeg[1]; // -1.14655

Expand All @@ -299,22 +302,24 @@ void Dune35tSimChannelExtractService::init() {
ycomb12 = ycomb11 + (ycomb11 - ycomb10);

// APA 2

geo->WireEndPoints(0,4,0,0,xyzbeg,xyzend); // first U wire in TPC 4.
constexpr geo::TPCID apa2{0, 4};
constexpr geo::PlaneID apa2_plane0{apa2,0};
geo->WireEndPoints(geo::WireID{apa2_plane0,0},xyzbeg,xyzend); // first U wire in TPC 4.
zcomb8 = xyzbeg[2]; // 51.924 -- same as above
ycomb17 = xyzend[1]; // 113.142 -- same as above

lastwire = geo->Nwires(0,4,0)-1; // 235 in v3
geo->WireEndPoints(0,4,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 4.
lastwire = geo->Nwires(apa2_plane0)-1; // 235 in v3
geo->WireEndPoints(geo::WireID{apa2_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 4.
zcomb11 = xyzend[2]; // 102.817 -- same as above
ycomb14 = xyzbeg[1]; // 0.83105

geo->WireEndPoints(0,4,1,0,xyzbeg,xyzend); // first V wire in TPC 4.
constexpr geo::PlaneID apa2_plane1{apa2, 1};
geo->WireEndPoints(geo::WireID{apa2_plane1,0},xyzbeg,xyzend); // first V wire in TPC 4.
zcomb10 = xyzend[2]; // 102.501 -- same as above
ycomb16 = xyzbeg[1]; // 113.142 -- everything ends here in y

lastwire = geo->Nwires(1,4,0)-1; // 227 in v3
geo->WireEndPoints(0,4,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 4.
lastwire = geo->Nwires(apa2_plane1)-1; // 227 in v3
geo->WireEndPoints(geo::WireID{apa2_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 4.
zcomb9 = xyzbeg[2]; // 52.2395 -- same as above
ycomb15 = xyzend[1]; // 1.14655

Expand All @@ -324,22 +329,24 @@ void Dune35tSimChannelExtractService::init() {
ycomb18 = ycomb17 + (ycomb17 - ycomb16);

// APA 3 -- a lot like APA 0

geo->WireEndPoints(0,6,0,0,xyzbeg,xyzend); // first U wire in TPC 6.
constexpr geo::TPCID apa3{0, 6};
constexpr geo::PlaneID apa3_plane0{apa3,0};
geo->WireEndPoints(geo::WireID{apa3_plane0,0},xyzbeg,xyzend); // first U wire in TPC 6.
zcomb14 = xyzbeg[2]; // 103.84
ycomb5 = xyzend[1]; // 113.142 -- same as above

lastwire = geo->Nwires(0,6,0)-1; // 358 in v3
geo->WireEndPoints(0,6,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 6.
lastwire = geo->Nwires(apa3_plane0)-1; // 358 in v3
geo->WireEndPoints(geo::WireID{apa3_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 6.
zcomb17 = xyzend[2]; // 154.741
ycomb2 = xyzbeg[1]; // -82.9389 -- same as above

geo->WireEndPoints(0,6,1,0,xyzbeg,xyzend); // first V wire in TPC 6.
constexpr geo::PlaneID apa3_plane1{apa3, 1};
geo->WireEndPoints(geo::WireID{apa3_plane1,0},xyzbeg,xyzend); // first V wire in TPC 6.
zcomb16 = xyzend[2]; // 154.425
ycomb4 = xyzbeg[1]; // 113.142 -- same as above

lastwire = geo->Nwires(1,6,0)-1; // 344 in v3
geo->WireEndPoints(0,6,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 6.
lastwire = geo->Nwires(apa3_plane1)-1; // 344 in v3
geo->WireEndPoints(geo::WireID{apa3_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 6.
zcomb15 = xyzbeg[2]; // 104.164
ycomb3 = xyzend[1]; // -82.6234 -- same as above

Expand Down
66 changes: 37 additions & 29 deletions dunesim/DetSim/SimWireDUNE35t_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -401,27 +401,29 @@ namespace detsim {
// wire running over the comb are beyond the endpoints. So we need to extrapolate.

double xyzbeg[3],xyzend[3];
int lastwire = 0;
unsigned int lastwire = 0;

// Numbers in comments are from Geometry V3 for debugging purposes.

// APA 0

geo->WireEndPoints(0,0,0,0,xyzbeg,xyzend); // first U wire in TPC 0.
constexpr geo::TPCID apa0{0,0};
constexpr geo::PlaneID apa0_plane0{apa0,0};
geo->WireEndPoints(geo::WireID{apa0_plane0,0},xyzbeg,xyzend); // first U wire in TPC 0.
zcomb2 = xyzbeg[2]; // 0.0
ycomb5 = xyzend[1]; // 113.142

lastwire = geo->Nwires(0,0,0)-1; // 358 in v3
geo->WireEndPoints(0,0,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 0.
lastwire = geo->Nwires(apa0_plane0)-1; // 358 in v3
geo->WireEndPoints(geo::WireID{apa0_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 0.
zcomb5 = xyzend[2]; // 50.8929
ycomb2 = xyzbeg[1]; // -82.9389

geo->WireEndPoints(0,0,1,0,xyzbeg,xyzend); // first V wire in TPC 0.
constexpr geo::PlaneID apa0_plane1{apa0, 1};
geo->WireEndPoints(geo::WireID{apa0_plane1,0},xyzbeg,xyzend); // first V wire in TPC 0.
zcomb4 = xyzend[2]; // 50.5774
ycomb4 = xyzbeg[1]; // 113.142

lastwire = geo->Nwires(1,0,0)-1; // 344 in v3
geo->WireEndPoints(0,0,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 0.
lastwire = geo->Nwires(apa0_plane1)-1; // 344 in v3
geo->WireEndPoints(geo::WireID{apa0_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 0.
zcomb3 = xyzbeg[2]; // 0.3155
ycomb3 = xyzend[1]; // -82.6234

Expand All @@ -439,22 +441,24 @@ namespace detsim {


// APA 1

geo->WireEndPoints(0,2,0,0,xyzbeg,xyzend); // first U wire in TPC 2.
constexpr geo::TPCID apa1{0, 2};
constexpr geo::PlaneID apa1_plane0{apa1,0};
geo->WireEndPoints(geo::WireID{apa1_plane0,0},xyzbeg,xyzend); // first U wire in TPC 2.
zcomb11 = xyzend[2]; // 102.817
ycomb8 = xyzbeg[1]; // -85.221

lastwire = geo->Nwires(0,2,0)-1; // 194 in v3
geo->WireEndPoints(0,2,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 2.
lastwire = geo->Nwires(apa1_plane0)-1; // 194 in v3
geo->WireEndPoints(geo::WireID{apa1_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 2.
zcomb8 = xyzbeg[2]; // 51.924
ycomb11 = xyzend[1]; // -0.831

geo->WireEndPoints(0,2,1,0,xyzbeg,xyzend); // first V wire in TPC 2.
constexpr geo::PlaneID apa1_plane1{apa1, 1};
geo->WireEndPoints(geo::WireID{apa1_plane1,0},xyzbeg,xyzend); // first V wire in TPC 2.
zcomb9 = xyzbeg[2]; // 52.2395
ycomb9 = xyzend[1]; // -85.222

lastwire = geo->Nwires(1,2,0)-1; // 188 in v3
geo->WireEndPoints(0,2,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 2.
lastwire = geo->Nwires(apa1_plane1)-1; // 188 in v3
geo->WireEndPoints(geo::WireID{apa1_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 2.
zcomb10 = xyzend[2]; // 102.501
ycomb10 = xyzbeg[1]; // -1.14655

Expand All @@ -470,22 +474,24 @@ namespace detsim {
ycomb12 = ycomb11 + (ycomb11 - ycomb10);

// APA 2

geo->WireEndPoints(0,4,0,0,xyzbeg,xyzend); // first U wire in TPC 4.
constexpr geo::TPCID apa2{0, 4};
constexpr geo::PlaneID apa2_plane0{apa2,0};
geo->WireEndPoints(geo::WireID{apa2_plane0,0},xyzbeg,xyzend); // first U wire in TPC 4.
zcomb8 = xyzbeg[2]; // 51.924 -- same as above
ycomb17 = xyzend[1]; // 113.142 -- same as above

lastwire = geo->Nwires(0,4,0)-1; // 235 in v3
geo->WireEndPoints(0,4,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 4.
lastwire = geo->Nwires(apa2_plane0)-1; // 235 in v3
geo->WireEndPoints(geo::WireID{apa2_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 4.
zcomb11 = xyzend[2]; // 102.817 -- same as above
ycomb14 = xyzbeg[1]; // 0.83105

geo->WireEndPoints(0,4,1,0,xyzbeg,xyzend); // first V wire in TPC 4.
constexpr geo::PlaneID apa2_plane1{apa2, 1};
geo->WireEndPoints(geo::WireID{apa2_plane1,0},xyzbeg,xyzend); // first V wire in TPC 4.
zcomb10 = xyzend[2]; // 102.501 -- same as above
ycomb16 = xyzbeg[1]; // 113.142 -- everything ends here in y

lastwire = geo->Nwires(1,4,0)-1; // 227 in v3
geo->WireEndPoints(0,4,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 4.
lastwire = geo->Nwires(apa2_plane1)-1; // 227 in v3
geo->WireEndPoints(geo::WireID{apa2_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 4.
zcomb9 = xyzbeg[2]; // 52.2395 -- same as above
ycomb15 = xyzend[1]; // 1.14655

Expand All @@ -502,22 +508,24 @@ namespace detsim {
ycomb18 = ycomb17 + (ycomb17 - ycomb16);

// APA 3 -- a lot like APA 0

geo->WireEndPoints(0,6,0,0,xyzbeg,xyzend); // first U wire in TPC 6.
constexpr geo::TPCID apa3{0, 6};
constexpr geo::PlaneID apa3_plane0{apa3,0};
geo->WireEndPoints(geo::WireID{apa3_plane0,0},xyzbeg,xyzend); // first U wire in TPC 6.
zcomb14 = xyzbeg[2]; // 103.84
ycomb5 = xyzend[1]; // 113.142 -- same as above

lastwire = geo->Nwires(0,6,0)-1; // 358 in v3
geo->WireEndPoints(0,6,0,lastwire,xyzbeg,xyzend); // last U wire in TPC 6.
lastwire = geo->Nwires(apa3_plane0)-1; // 358 in v3
geo->WireEndPoints(geo::WireID{apa3_plane0,lastwire},xyzbeg,xyzend); // last U wire in TPC 6.
zcomb17 = xyzend[2]; // 154.741
ycomb2 = xyzbeg[1]; // -82.9389 -- same as above

geo->WireEndPoints(0,6,1,0,xyzbeg,xyzend); // first V wire in TPC 6.
constexpr geo::PlaneID apa3_plane1{apa3, 1};
geo->WireEndPoints(geo::WireID{apa3_plane1,0},xyzbeg,xyzend); // first V wire in TPC 6.
zcomb16 = xyzend[2]; // 154.425
ycomb4 = xyzbeg[1]; // 113.142 -- same as above

lastwire = geo->Nwires(1,6,0)-1; // 344 in v3
geo->WireEndPoints(0,6,1,lastwire,xyzbeg,xyzend); // last V wire in TPC 6.
lastwire = geo->Nwires(apa3_plane1)-1; // 344 in v3
geo->WireEndPoints(geo::WireID{apa3_plane1,lastwire},xyzbeg,xyzend); // last V wire in TPC 6.
zcomb15 = xyzbeg[2]; // 104.164
ycomb3 = xyzend[1]; // -82.6234 -- same as above

Expand Down
28 changes: 13 additions & 15 deletions dunesim/DetectorVariations/RandomChannelStatusService_service.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -73,15 +73,15 @@ namespace detvar
{
public:
RandomTPC(const geo::GeometryCore* geom)
: fCryos(VectorViaSet<geo::CryostatID>(geom->begin_cryostat_id(),
geom->end_cryostat_id()))
: fCryos(VectorViaSet<geo::CryostatID>(geom->begin<geo::CryostatID>(),
geom->end<geo::CryostatID>()))
{
for(geo::CryostatID c: fCryos){
fTPCs[c] = VectorViaSet<geo::TPCID>(geom->begin_TPC_id(c),
geom->end_TPC_id(c));
fTPCs[c] = VectorViaSet<geo::TPCID>(geom->begin<geo::TPCID>(c),
geom->end<geo::TPCID>(c));

fTPCsets[c] = VectorViaSet<readout::TPCsetID>(geom->begin_TPCset_id(c),
geom->end_TPCset_id(c));
fTPCsets[c] = VectorViaSet<readout::TPCsetID>(geom->begin<readout::TPCsetID>(c),
geom->end<readout::TPCsetID>(c));
}
}

Expand Down Expand Up @@ -149,7 +149,7 @@ namespace detvar
// No good way of enumerating all the unique channels, or selecting them by
// index. Have to figure them out from the wires.
std::set<raw::ChannelID_t> chanset[3];
for(geo::WireID wire: geom->IterateWireIDs(tpc)){
for(geo::WireID const& wire: geom->Iterate<geo::WireID>(tpc)){
const raw::ChannelID_t chan = geom->PlaneWireToChannel(wire);
// But this also gives us wires that are actually attached to the other
// face and just wrapped onto this face. So long as the order of the
Expand Down Expand Up @@ -217,7 +217,7 @@ namespace detvar

// goodchans = allchans - badchans
std::set<raw::ChannelID_t> allchans;
for(geo::WireID wire: geom->IterateWireIDs())
for(geo::WireID const& wire: geom->Iterate<geo::WireID>())
allchans.insert(geom->PlaneWireToChannel(wire));

std::set_difference(allchans.begin(), allchans.end(),
Expand All @@ -233,7 +233,7 @@ namespace detvar
// Geometry doesn't have a way to iterate directly over channels. Iterate
// over the wires and convert them. Use a set to remove duplicates
std::set<raw::ChannelID_t> allchans;
for(geo::WireID wire: geom->IterateWireIDs())
for(geo::WireID const& wire: geom->Iterate<geo::WireID>())
allchans.insert(geom->PlaneWireToChannel(wire));

// But a vector is much easier to pick from randomly. This will be used for
Expand All @@ -258,11 +258,11 @@ namespace detvar
art::ServiceHandle<geo::Geometry> geom;

std::map<readout::TPCsetID, std::set<raw::ChannelID_t>> tpcset_to_chans;
for(const readout::TPCsetID& ts: geom->IterateTPCsetIDs()){
for(const readout::TPCsetID& ts: geom->Iterate<readout::TPCsetID>()){
// There is no version of IterateWireIDs over a TPCset. Use another layer
// of indirection.
for(geo::TPCID t: geom->TPCsetToTPCs(ts)){
for(const geo::WireID& wire: geom->IterateWireIDs(t)){
for(const geo::WireID& wire: geom->Iterate<geo::WireID>(t)){
const raw::ChannelID_t chan = geom->PlaneWireToChannel(wire);
tpcset_to_chans[ts].insert(chan);
}
Expand Down Expand Up @@ -290,18 +290,16 @@ namespace detvar
art::ServiceHandle<geo::Geometry> geom;

std::map<geo::TPCID, std::vector<raw::ChannelID_t>> tpc_to_chans;
for(const geo::TPCID& t: geom->IterateTPCIDs()){
for(const geo::WireID& wire: geom->Iterate<geo::WireID>()){
// Geometry doesn't provide a way to directly iterate the channels in the
// TPC. Instead iterate the wires and convert to channels
for(const geo::WireID& wire: geom->IterateWireIDs(t)){
const raw::ChannelID_t chan = geom->PlaneWireToChannel(wire);
// But this also gives us wires that are actually attached to the other
// face and just wrapped onto this face. So long as the order of the
// vector returned from this function is meaningful, this should work
// to keep just the ones we need.
if(geom->ChannelToWire(chan)[0] == wire)
tpc_to_chans[t].push_back(chan);
}
tpc_to_chans[wire.asPlaneID().asTPCID()].push_back(chan);
}

RandomTPC tpcs(geom.get());
Expand Down
8 changes: 4 additions & 4 deletions dunesim/EventGenerator/CORSIKAprotodunedp/CORSIKAGendp_module.cc
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -373,9 +373,9 @@ namespace evgendp{

//compute shower area based on the maximal x,z dimensions of cryostat boundaries + fShowerAreaExtension
art::ServiceHandle<geo::Geometry> geom;
for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
for(auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
double bounds[6] = {0.};
geom->Cryostat(c).Boundaries(bounds);
cryostat.Boundaries(bounds);
for (unsigned int bnd = 0; bnd<6; bnd++){
mf::LogVerbatim("CORSIKAGendp")<<"Cryo Boundary: "<<bnd<<"="<<bounds[bnd]<<" ( + Buffer="<<fBuffBox[bnd]<<")\n";
if(fabs(bounds[bnd])>fabs(fShowerBounds[bnd])){
Expand Down Expand Up @@ -662,9 +662,9 @@ namespace evgendp{

// now check if the particle goes through any cryostat in the detector
// if so, add it to the truth object.
for(unsigned int c = 0; c < geom->Ncryostats(); ++c){
for(auto const& cryostat : geom->Iterate<geo::CryostatGeo>()) {
double bounds[6] = {0.};
geom->Cryostat(c).Boundaries(bounds);
cryostat.Boundaries(bounds);

//add a buffer box around the cryostat bounds to increase the acceptance and account for scattering
//By default, the buffer box has zero size
Expand Down