Skip to content
Closed
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 5 additions & 4 deletions sbncode/CAFMaker/CAFMaker_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -92,6 +92,7 @@
#include "lardataobj/RecoBase/PFParticle.h"
#include "lardataobj/RecoBase/Slice.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/TrackHitMeta.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/Shower.h"
#include "lardataobj/RecoBase/MCSFitResult.h"
Expand Down Expand Up @@ -1384,7 +1385,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
}

// Prep truth-to-reco-matching info
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map;
std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
std::map<int, caf::HitsEnergy> id_to_hit_energy_map;

Expand Down Expand Up @@ -1916,8 +1917,8 @@ void CAFMaker::produce(art::Event& evt) noexcept {
FindManyPStrict<recob::Vertex>(fmPFPart, evt,
fParams.PFParticleLabel() + slice_tag_suff);

art::FindManyP<recob::Hit> fmTrackHit =
FindManyPStrict<recob::Hit>(slcTracks, evt,
art::FindManyP<recob::Hit, recob::TrackHitMeta> fmTrackHit =
FindManyPDStrict<recob::Hit, recob::TrackHitMeta>(slcTracks, evt,
fParams.RecoTrackLabel() + slice_tag_suff);

art::FindManyP<recob::Hit> fmShowerHit =
Expand Down Expand Up @@ -2162,7 +2163,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
FillTrackDazzle(fmTrackDazzle.at(iPart).front(), trk);
}
if (fmCalo.isValid()) {
FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart),
FillTrackCalo(fmCalo.at(iPart), *thisTrack[0], fmTrackHit.at(iPart), fmTrackHit.data(iPart),
(fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(),
fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(),
dprop, trk);
Expand Down
47 changes: 45 additions & 2 deletions sbncode/CAFMaker/FillReco.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -10,6 +10,10 @@
#include "FillReco.h"
#include "RecoUtils/RecoUtils.h"

#include "larevt/SpaceCharge/SpaceCharge.h"
#include "larevt/SpaceChargeServices/SpaceChargeService.h"
#include "larcore/CoreUtils/ServiceUtil.h"

namespace caf
{

Expand Down Expand Up @@ -720,8 +724,28 @@ namespace caf
}
}

// Helper function: get the e field
double GetEfield(const detinfo::DetectorPropertiesData& dprop, const geo::Point_t loc) {
auto const* sce = lar::providerFrom<spacecharge::SpaceChargeService>();

double EField = dprop.Efield();
if (sce->EnableSimEfieldSCE()) {
// Gets relative E field Distortions
geo::Vector_t EFieldOffsets = sce->GetEfieldOffsets(loc);
// Add 1 in X direction as this is the direction of the drift field
EFieldOffsets = EFieldOffsets + geo::Vector_t{1, 0, 0};
// Convert to Absolute E Field from relative
EFieldOffsets = EField * EFieldOffsets;
// We only care about the magnitude for recombination
EField = EFieldOffsets.r();
}
return EField;
}

void FillTrackPlaneCalo(const anab::Calorimetry &calo,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrackCalo &srcalo) {
Expand Down Expand Up @@ -758,7 +782,8 @@ namespace caf

// lookup the wire -- the Calorimery object makes this
// __way__ harder than it should be
for (const art::Ptr<recob::Hit> &h: hits) {
for (unsigned i_hit = 0; i_hit < hits.size(); i_hit++) {
const art::Ptr<recob::Hit> &h = hits[i_hit];
if (h.key() == tps[i]) {
p.wire = h->WireID().Wire;
p.tpc = h->WireID().TPC;
Expand All @@ -770,6 +795,22 @@ namespace caf
p.mult = h->Multiplicity();
p.start = h->StartTick();
p.end = h->EndTick();


// Get the trajectory point index from this hit. Again -- this is too hard.
//
// Use this to get the (SCE corrected) efield and the angle to the drift direction
unsigned traj_point_index = thms.at(i_hit)->Index();
unsigned int int_max_as_unsigned_int{std::numeric_limits<int>::max()};
if (traj_point_index != int_max_as_unsigned_int && // invalid
track.HasValidPoint(traj_point_index)) {
float costh_drift = track.DirectionAtPoint(traj_point_index).X();
// p.costh_drift = costh_drift;
float efield = GetEfield(dprop, track.LocationAtPoint(traj_point_index));
// p.efield = efield;
(void) costh_drift;
(void) efield;
}
}
}

Expand Down Expand Up @@ -823,7 +864,9 @@ namespace caf
}

void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrack& srtrack,
Expand All @@ -838,7 +881,7 @@ namespace caf
if (calo.PlaneID()) {
unsigned plane_id = calo.PlaneID().Plane;
assert(plane_id < 3);
FillTrackPlaneCalo(calo, hits, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]);
FillTrackPlaneCalo(calo, track, hits, thms, fill_calo_points, fillhit_rrstart, fillhit_rrend, dprop, srtrack.calo[plane_id]);
}
}

Expand Down
5 changes: 5 additions & 0 deletions sbncode/CAFMaker/FillReco.h
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@
#include "lardataobj/RecoBase/Shower.h"
#include "lardataobj/RecoBase/Slice.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/TrackHitMeta.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/SpacePoint.h"
Expand Down Expand Up @@ -168,7 +169,9 @@ namespace caf
bool allowEmpty = false);

void FillTrackPlaneCalo(const anab::Calorimetry &calo,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrackCalo &srcalo);
Expand All @@ -186,7 +189,9 @@ namespace caf
bool allowEmpty = false);

void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
const recob::Track& track,
const std::vector<art::Ptr<recob::Hit>> &hits,
const std::vector<const recob::TrackHitMeta*>& thms,
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
const detinfo::DetectorPropertiesData &dprop,
caf::SRTrack& srtrack,
Expand Down
26 changes: 13 additions & 13 deletions sbncode/CAFMaker/FillTrue.cxx
Original file line number Diff line number Diff line change
Expand Up @@ -136,7 +136,7 @@ namespace caf {
}//FillTrackTruth

// Assumes truth matching and calo-points are filled
void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
void FillTrackCaloTruth(const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::vector<simb::MCParticle> &mc_particles,
const geo::GeometryCore& geometry,
const geo::WireReadoutGeom& wireReadout,
Expand All @@ -163,10 +163,10 @@ namespace caf {

// Load the hits
// match on the channel, which is unique
const std::vector<std::pair<geo::WireID, const sim::IDE*>> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
const std::vector<caf::ParticleIDE> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
std::map<unsigned, std::vector<const sim::IDE *>> chan_2_ides;
for (auto const &ide_pair: match_ides) {
chan_2_ides[wireReadout.PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second);
for (auto const &ide_p: match_ides) {
chan_2_ides[wireReadout.PlaneWireToChannel(ide_p.wire)].push_back(ide_p.ide);
}

// pre-compute partial ranges
Expand Down Expand Up @@ -634,15 +634,15 @@ namespace caf {
void FillTrueG4Particle(const simb::MCParticle &particle,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const cheat::BackTrackerService &backtracker,
const cheat::ParticleInventoryService &inventory_service,
const std::vector<art::Ptr<simb::MCTruth>> &neutrinos,
caf::SRTrueParticle &srparticle) {

std::vector<std::pair<geo::WireID, const sim::IDE *>> empty;
const std::vector<std::pair<geo::WireID, const sim::IDE *>> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;
std::vector<caf::ParticleIDE> empty;
const std::vector<caf::ParticleIDE> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;

std::vector<art::Ptr<recob::Hit>> emptyHits;
const std::vector<art::Ptr<recob::Hit>> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits;
Expand All @@ -662,9 +662,9 @@ namespace caf {
}
}

for (auto const &ide_pair: particle_ides) {
const geo::WireID &w = ide_pair.first;
const sim::IDE *ide = ide_pair.second;
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 && w.Plane < 3 && w.Cryostat < 2){
srparticle.plane[w.Cryostat][w.Plane].visE += ide->energy / 1000. /* MeV -> GeV*/;
Expand Down Expand Up @@ -882,8 +882,8 @@ namespace caf {
return ret;
}

std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> ret;
std::map<int, std::vector<caf::ParticleIDE>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout) {
std::map<int, std::vector<caf::ParticleIDE>> ret;

for (const art::Ptr<sim::SimChannel> sc : simchannels) {
// Lookup the wire of this channel
Expand All @@ -895,7 +895,7 @@ namespace caf {
for (const auto &item : sc->TDCIDEMap()) {
for (const sim::IDE &ide: item.second) {
// indexing initializes empty vector
ret[abs(ide.trackID)].push_back({thisWire, &ide});
ret[abs(ide.trackID)].push_back({thisWire, item.first, &ide});
}
}
}
Expand Down
13 changes: 10 additions & 3 deletions sbncode/CAFMaker/FillTrue.h
Original file line number Diff line number Diff line change
Expand Up @@ -40,6 +40,13 @@ namespace caf
float totE;
};

struct ParticleIDE {
geo::WireID wire;
unsigned short tick;
const sim::IDE *ide;
};


// Helpers
caf::Wall_t GetWallCross( const geo::BoxBoundedGeo &volume,
const TVector3 p0,
Expand Down Expand Up @@ -73,7 +80,7 @@ namespace caf
void FillTrueG4Particle(const simb::MCParticle &particle,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const cheat::BackTrackerService &backtracker,
const cheat::ParticleInventoryService &inventory_service,
Expand Down Expand Up @@ -103,7 +110,7 @@ namespace caf
caf::SRTrack& srtrack,
bool allowEmpty = false);

void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
void FillTrackCaloTruth(const std::map<int, std::vector<ParticleIDE>> &id_to_ide_map,
const std::vector<simb::MCParticle> &mc_particles,
const geo::GeometryCore & geometry,
const geo::WireReadoutGeom& wireReadout,
Expand Down Expand Up @@ -133,7 +140,7 @@ namespace caf
CLHEP::HepRandomEngine &rand,
std::vector<caf::SRFakeReco> &srfakereco);

std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout);
std::map<int, std::vector<ParticleIDE>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::WireReadoutGeom &wireReadout);
std::map<int, std::vector<art::Ptr<recob::Hit>>> PrepTrueHits(const std::vector<art::Ptr<recob::Hit>> &allHits,
const detinfo::DetectorClocksData &clockData, const cheat::BackTrackerService &backtracker);
std::map<int, caf::HitsEnergy> SetupIDHitEnergyMap(const std::vector<art::Ptr<recob::Hit>> &allHits, const detinfo::DetectorClocksData &clockData,
Expand Down
3 changes: 2 additions & 1 deletion sbncode/Calibration/TrackCaloSkimmer.h
Original file line number Diff line number Diff line change
Expand Up @@ -64,6 +64,7 @@
#include "sbnobj/Common/CRT/CRTHitT0TaggingInfo.cc"
#include "sbnobj/Common/CRT/CRTHitT0TaggingTruthInfo.cc"

#include "sbncode/CAFMaker/FillTrue.h"
#include "ITCSSelectionTool.h"

namespace sbn {
Expand Down Expand Up @@ -161,7 +162,7 @@ class sbn::TrackCaloSkimmer : public art::EDAnalyzer {
const std::vector<art::Ptr<simb::MCParticle>> &mcparticles,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
Expand Down
33 changes: 17 additions & 16 deletions sbncode/Calibration/TrackCaloSkimmer_module.cc
Original file line number Diff line number Diff line change
Expand Up @@ -289,7 +289,7 @@ void sbn::TrackCaloSkimmer::analyze(art::Event const& e)
// Prep truth-to-reco-matching info
//
// Use helper functions from CAFMaker/FillTrue
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map;
std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map;
std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map;
const cheat::BackTrackerService *bt = NULL;

Expand Down Expand Up @@ -545,14 +545,14 @@ geo::Point_t WireToTrajectoryPosition(const geo::Point_t &loc, const geo::TPCID
sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE *>>> &id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> &id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> &id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
const geo::WireReadoutGeom *wireReadout) {

std::vector<std::pair<geo::WireID, const sim::IDE *>> empty;
const std::vector<std::pair<geo::WireID, const sim::IDE *>> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;
std::vector<caf::ParticleIDE> empty;
const std::vector<caf::ParticleIDE> &particle_ides = id_to_ide_map.count(particle.TrackId()) ? id_to_ide_map.at(particle.TrackId()) : empty;

std::vector<art::Ptr<recob::Hit>> emptyHits;
const std::vector<art::Ptr<recob::Hit>> &particle_hits = id_to_truehit_map.count(particle.TrackId()) ? id_to_truehit_map.at(particle.TrackId()) : emptyHits;
Expand All @@ -569,9 +569,9 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
trueparticle.plane0nhit = 0;
trueparticle.plane1nhit = 0;
trueparticle.plane2nhit = 0;
for (auto const &ide_pair: particle_ides) {
const geo::WireID &w = ide_pair.first;
const sim::IDE *ide = ide_pair.second;
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*/;
Expand Down Expand Up @@ -720,10 +720,11 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
// Organize deposition info into per-wire true "Hits" -- key is the Channel Number
std::map<unsigned, sbn::TrueHit> truehits;

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;
unsigned c = wireReadout->PlaneWireToChannel(w);
const sim::IDE *ide = ide_pair.second;
const sim::IDE *ide = ide_part.ide;
unsigned short tick = ide_part.tick;

// Set stuff
truehits[c].cryo = w.Cryostat;
Expand All @@ -739,6 +740,8 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
truehits[c].p.y = (truehits[c].p.y*old_elec + ide->y*ide->numElectrons) / new_elec;
truehits[c].p.z = (truehits[c].p.z*old_elec + ide->z*ide->numElectrons) / new_elec;

truehits[c].time = (truehits[c].time*old_elec + tick*ide->numElectrons) / new_elec;

// Also get the position with space charge un-done
geo::Point_t ide_p(ide->x, ide->y, ide->z);
geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w);
Expand All @@ -754,10 +757,10 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,
}

// Compute widths
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;
unsigned c = wireReadout->PlaneWireToChannel(w);
const sim::IDE *ide = ide_pair.second;
const sim::IDE *ide = ide_part.ide;

geo::Point_t ide_p(ide->x, ide->y, ide->z);
geo::Point_t ide_p_scecorr = WireToTrajectoryPosition(ide_p, w);
Expand All @@ -782,8 +785,6 @@ sbn::TrueParticle TrueParticleInfo(const simb::MCParticle &particle,

// Compute the time of each hit
for (sbn::TrueHit &h: truehits_v) {
h.time = dprop.ConvertXToTicks(h.p.x, h.plane, h.tpc, h.cryo);

double xdrift = abs(h.p.x - wireReadout->Plane(geo::PlaneID(h.cryo, h.tpc, 0)).GetCenter().X());
h.tdrift = xdrift / dprop.DriftVelocity();
}
Expand Down Expand Up @@ -1010,7 +1011,7 @@ void sbn::TrackCaloSkimmer::FillTrackTruth(const detinfo::DetectorClocksData &cl
const std::vector<art::Ptr<simb::MCParticle>> &mcparticles,
const std::vector<geo::BoxBoundedGeo> &active_volumes,
const std::vector<std::vector<geo::BoxBoundedGeo>> &tpc_volumes,
const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> id_to_ide_map,
const std::map<int, std::vector<caf::ParticleIDE>> id_to_ide_map,
const std::map<int, std::vector<art::Ptr<recob::Hit>>> id_to_truehit_map,
const detinfo::DetectorPropertiesData &dprop,
const geo::GeometryCore *geo,
Expand Down
1 change: 1 addition & 0 deletions sbncode/EventGenerator/CMakeLists.txt
Original file line number Diff line number Diff line change
@@ -1,2 +1,3 @@
add_subdirectory(MeVPrtl)
add_subdirectory(MultiPart)
add_subdirectory(WireModGen)
Loading