WireMod for Develop#650
Conversation
[Bugfix] NC Delta Rad Filter was not working. Copied working RES structure and changed to filter to NC Delta Rad
… into production/sbnd-gen2
gputnam
left a comment
There was a problem hiding this comment.
Looks good! Insofar as it intersects with stuff I've contributed. It looks like a SinglePhoton analysis commit got in here somehow though? Bit confused by that. But I don't see that as holding this up.
|
The WireMod utility portions look good to me. However, I do still have some concerns about the |
|
Converted to a draft while waiting for #654 to be merged. |
PetrilloAtWork
left a comment
There was a problem hiding this comment.
Since @hausnerh asked for a comment, I have decorated the module with a quick analysis of the issue.
The question is whether any ambiguity from ChannelToWire() result might create trouble.
My findings:
- there is the potential for ambiguity, since that function can and will return multiple wires. The wires will all belong to the same cryostat and plane, but will have different wire number and may have (usually do have) different TPC number;
- the code uses
sbn::ReadoutIDEto carry around the information, and that one does not have room for a channel number, only fora wire ID; - most of the occurrences of the use of this wire ID are to select the plane number, which is expected to be always correct;
- one exception is the use of
WireToTrajectoryPosition()which, despite the name, just corrects a 3D point for field distortions; this function uses only the TPC ID, not the wire number. If I recall correctly, the TPC ID is used only to determine the drift direction (after all, we know the exact location of the point). While the TPC ID may be wrong, in practice in SBND all TPCs sharing a channel have the same drift direction, so this does not become an issue (DUNE APA would be a different story). - the other exception is the wire and TPC number values being stored in the output. This is the only instance where the effects of the potentially wrong information will be observable. I have left a possible workaround.
Do with it what you will :-D
| const geo::WireID &w = ide_p.wire; | ||
| const sim::IDE *ide = ide_p.ide; | ||
|
|
||
| if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2){ | ||
| srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; |
There was a problem hiding this comment.
To help addressing @hausnerh's concern: here the wire number, which is the only ambiguous part when going from a channel to a wire ID, is not used. I suggest this code makes it explicit:
| const geo::WireID &w = ide_p.wire; | |
| const sim::IDE *ide = ide_p.ide; | |
| if(w.Plane >= 0 && w.Plane < 3 && w.Cryostat < 2){ | |
| srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; | |
| const geo::PlaneID &p = ide_p.wire.asPlaneID(); | |
| const sim::IDE *ide = ide_p.ide; | |
| if(p.Plane >= 0 && p.Plane < 3 && p.Cryostat < 2){ | |
| srparticle.plane[p.Cryostat][p.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/; |
(.asPlaneID() is not necessary since conversion can be implicit, but it helps the reader)
| for (auto const &ide_p: particle_ides) { | ||
| const geo::WireID &w = ide_p.wire; | ||
| const sim::IDE *ide = ide_p.ide; | ||
|
|
||
| if (w.Plane == 0) { | ||
| trueparticle.plane0VisE += ide->energy / 1000. /* MeV -> GeV*/; | ||
| } | ||
| else if (w.Plane == 1) { | ||
| trueparticle.plane1VisE += ide->energy / 1000. /* MeV -> GeV*/; | ||
| } | ||
| else if (w.Plane == 2) { | ||
| trueparticle.plane2VisE += ide->energy / 1000. /* MeV -> GeV*/; | ||
| } | ||
| } |
There was a problem hiding this comment.
Here only the plane number is used:
| for (auto const &ide_p: particle_ides) { | |
| const geo::PlaneID::Plane_t p = ide_p.wire.Plane; | |
| const sim::IDE *ide = ide_p.ide; | |
| switch (p) { | |
| case 0: | |
| trueparticle.plane0VisE += ide->energy / 1000. /* MeV -> GeV*/; | |
| break; | |
| case 1: | |
| trueparticle.plane1VisE += ide->energy / 1000. /* MeV -> GeV*/; | |
| break; | |
| case 2: | |
| trueparticle.plane2VisE += ide->energy / 1000. /* MeV -> GeV*/; | |
| break; | |
| } | |
| } |
If you don't care of that level of pedantry, an int instead of the format geo::PlaneID::Plane_t will do just as well.
I find a switch more appropriate here, but that's not really the point of this comment.
| for (auto const &ide_part: particle_ides) { | ||
| const geo::WireID &w = ide_part.wire; |
There was a problem hiding this comment.
Here:
PlaneWireToChannel()will correctly track back to the channel, no matter which of its wire was chosen, so it's not a problem;WireToTrajectoryPosition()(whose name is wrong) uses only the TPC information, so one can make that explicit too:
| for (auto const &ide_part: particle_ides) { | |
| const geo::WireID &w = ide_part.wire; | |
| for (auto const &ide_part: particle_ides) { | |
| const geo::WireID &w = ide_part.wire; | |
| const geo::TPCID &t = w.asTPCID(); |
and then use t in the following call.
| for (auto const &ide_pair: particle_ides) { | ||
| const geo::WireID &w = ide_pair.first; | ||
| for (auto const &ide_part: particle_ides) { | ||
| const geo::WireID &w = ide_part.wire; |
There was a problem hiding this comment.
This can be treated as the case below (sorry for the back-and-forth)...
| const geo::WireID &w = ide_part.wire; | |
| const geo::WireID &w = ide_part.wire; | |
| const geo::TPCID &t = w.asTPCID(); |
with one exception.
| truehits[c].cryo = w.Cryostat; | ||
| truehits[c].tpc = w.TPC; | ||
| truehits[c].plane = w.Plane; | ||
| truehits[c].wire = w.Wire; |
There was a problem hiding this comment.
This one looks like the only case where the wrong wire number may have an effect.
If this is unacceptable, then one (expensive) workaround is to extracted the wire again:
// track down the actual TPC and wire for the IDE
geo::Point_t const IDEpos{ ide->x, ide->y, ide->z };
geo::TPCID const IDETPC = wireReadout->FindTPCatPosition(IDEpos);
geo::WireID const IDEwire = IDETPC.isValid
? wireReadout->NearestWireID(IDEpos, { IDETPC, w.Plane }): geo::WireID{};Note that we still select the relevant plane number from the cached wire.
Do this above, and then you can use IDETPC for every time you need a TPC ID (if you need a plane ID too, better to define it with geo::planeID const IDEplane{ IDETPC, w.Plane };).
However, note that this returns the real TPC and wire of the IDE; if that IDE was shifted to simulate some displacement, its wire may end up being a different one than w, and on a different channel than c.
| raw::ChannelID_t channel = sc->Channel(); | ||
| std::vector<geo::WireID> maybewire = wireReadout.ChannelToWire(channel); | ||
| geo::WireID thisWire; // Default constructor makes invalid wire | ||
| if (maybewire.size()) thisWire = maybewire[0]; |
There was a problem hiding this comment.
Considering that
- often this value is not checked downstream
- a
SimChannelshould not be associated to wireless channels
I suggest to make this function loudly complain in case of inconsistency:
| if (maybewire.size()) thisWire = maybewire[0]; | |
| if (maybewire.empty()) { | |
| if (sc->TDCIDEMap().empty()) continue; | |
| throw cet::exception("sbn::PrepSimChannels") | |
| << "Data on sim::SimChannel CH=" << channel | |
| << ", which has no associated wires!\n"; | |
| } | |
| thisWire = maybewire[0]; |
Description
Add WireMod utility functions to
develop. This also brings in refactor of RecoUtils. The WireMod utility should be identical to the SBND Gen 1 production based on v10_06_00_10. This is the corresponding PR to sbncode #649.the description.