Skip to content

Commit de192fe

Browse files
ibsafaibsafa
authored
Merge branch 'develop' into feature/SBNDBNB_Retriever
2 parents 83f8895 + d6f3a30 commit de192fe

Some content is hidden

Large Commits have some content hidden by default. Use the searchbox below for content that may be hidden.

59 files changed

+631
-452
lines changed

‎CMakeLists.txt‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -16,7 +16,7 @@
1616
cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
1717

1818
find_package(cetmodules 3.20.00 REQUIRED)
19-
project(sbncode VERSION 09.91.02.01 LANGUAGES CXX)
19+
project(sbncode VERSION 10.01.03 LANGUAGES CXX)
2020

2121
message(STATUS "\n\n ========================== ${PROJECT_NAME} ==========================")
2222

‎sbncode/CAFMaker/CAFMaker_module.cc‎

Lines changed: 16 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -87,6 +87,8 @@
8787
#include "messagefacility/MessageLogger/MessageLogger.h"
8888

8989
// LArSoft includes
90+
#include "larcore/Geometry/WireReadout.h"
91+
#include "larcore/Geometry/Geometry.h"
9092
#include "lardataobj/RecoBase/PFParticle.h"
9193
#include "lardataobj/RecoBase/Slice.h"
9294
#include "lardataobj/RecoBase/Track.h"
@@ -99,8 +101,6 @@
99101
#include "nusimdata/SimulationBase/MCNeutrino.h"
100102
#include "nusimdata/SimulationBase/GTruth.h"
101103

102-
#include "fhiclcpp/ParameterSetRegistry.h"
103-
104104
#include "sbnobj/Common/EventGen/MeVPrtl/MeVPrtlTruth.h"
105105
#include "sbnobj/Common/Reco/RangeP.h"
106106
#include "sbnobj/Common/SBNEventWeight/EventWeightMap.h"
@@ -1302,7 +1302,9 @@ void CAFMaker::produce(art::Event& evt) noexcept {
13021302
auto const clock_data = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
13031303
auto const dprop =
13041304
art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clock_data);
1305-
const geo::GeometryCore *geometry = lar::providerFrom<geo::Geometry>();
1305+
const geo::GeometryCore* geom = lar::providerFrom<geo::Geometry>();
1306+
const geo::WireReadoutGeom &wireReadout =
1307+
art::ServiceHandle<geo::WireReadout>()->Get();
13061308

13071309
auto const *sce = lar::providerFrom<spacecharge::SpaceChargeService>();
13081310

@@ -1330,7 +1332,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
13301332
if ( !isRealData ) {
13311333
art::ServiceHandle<cheat::BackTrackerService> bt_serv;
13321334

1333-
id_to_ide_map = PrepSimChannels(simchannels, *geometry);
1335+
id_to_ide_map = PrepSimChannels(simchannels, wireReadout);
13341336
id_to_truehit_map = PrepTrueHits(hits, clock_data, *bt_serv);
13351337
id_to_hit_energy_map = SetupIDHitEnergyMap(hits, clock_data, *bt_serv);
13361338
}
@@ -1550,22 +1552,21 @@ void CAFMaker::produce(art::Event& evt) noexcept {
15501552
}
15511553
}
15521554

1553-
// Get all of the CRTPMT Matches ..
1555+
// Get all of the CRTPMT Matches
15541556
std::vector<caf::SRCRTPMTMatch> srcrtpmtmatches;
1555-
std::cout << "srcrtpmtmatches.size = " << srcrtpmtmatches.size() << "\n";
15561557
art::Handle<std::vector<sbn::crt::CRTPMTMatching>> crtpmtmatch_handle;
15571558
GetByLabelStrict(evt, fParams.CRTPMTLabel(), crtpmtmatch_handle);
15581559
if(crtpmtmatch_handle.isValid()){
1559-
std::cout << "valid handle! label: " << fParams.CRTPMTLabel() << "\n";
15601560
const std::vector<sbn::crt::CRTPMTMatching> &crtpmtmatches = *crtpmtmatch_handle;
15611561
for (unsigned i = 0; i < crtpmtmatches.size(); i++) {
1562+
int topen = 0, topex = 0, sideen = 0, sideex = 0; // bottomen = 0, bottomex = 0;
15621563
srcrtpmtmatches.emplace_back();
1563-
FillCRTPMTMatch(crtpmtmatches[i],srcrtpmtmatches.back());
1564+
FillCRTPMTMatch(crtpmtmatches[i],
1565+
topen, topex, sideen, sideex,
1566+
//bottomen, bottomex,
1567+
srcrtpmtmatches.back());
15641568
}
15651569
}
1566-
else{
1567-
std::cout << "crtpmtmatch_handle.isNOTValid!\n";
1568-
}
15691570

15701571
// Get all of the OpFlashes
15711572
std::vector<caf::SROpFlash> srflashes;
@@ -2048,7 +2049,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
20482049
FillTrackRangeP(*thisTrack[0], rangePs, trk);
20492050

20502051
if (fmChi2PID.isValid()) {
2051-
FillTrackChi2PID(fmChi2PID.at(iPart), lar::providerFrom<geo::Geometry>(), trk);
2052+
FillTrackChi2PID(fmChi2PID.at(iPart), trk);
20522053
}
20532054
if (fmScatterClosestApproach.isValid() && fmScatterClosestApproach.at(iPart).size()==1) {
20542055
FillTrackScatterClosestApproach(fmScatterClosestApproach.at(iPart).front(), trk);
@@ -2063,7 +2064,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
20632064
FillTrackCalo(fmCalo.at(iPart), fmTrackHit.at(iPart),
20642065
(fParams.FillHitsNeutrinoSlices() && NeutrinoSlice) || fParams.FillHitsAllSlices(),
20652066
fParams.TrackHitFillRRStartCut(), fParams.TrackHitFillRREndCut(),
2066-
lar::providerFrom<geo::Geometry>(), dprop, trk);
2067+
dprop, trk);
20672068
}
20682069
if (fmCRTHitMatch.isValid() && fDet == kICARUS) {
20692070
art::FindManyP<sbn::crt::CRTHit> CRTT02Hit = FindManyPStrict<sbn::crt::CRTHit>
@@ -2101,7 +2102,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
21012102
FillTrackTruth(fmTrackHit.at(iPart), id_to_hit_energy_map, true_particles, clock_data, trk);
21022103
// Hit truth information corresponding to Calo-Points
21032104
// Assumes truth matching and calo-points are filled
2104-
if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, geometry, clock_data, sce, trk);
2105+
if (mc_particles.isValid() && fParams.FillTrackCaloTruth()) FillTrackCaloTruth(id_to_ide_map, *mc_particles, *geom, wireReadout, clock_data, sce, trk);
21052106
}
21062107
}
21072108
} // thisTrack exists
@@ -2110,7 +2111,7 @@ void CAFMaker::produce(art::Event& evt) noexcept {
21102111
assert(thisShower.size() == 1);
21112112

21122113
SRShower& shw = pfp.shw;
2113-
FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), lar::providerFrom<geo::Geometry>(), producer, shw);
2114+
FillShowerVars(*thisShower[0], vertex, fmShowerHit.at(iPart), wireReadout, producer, shw);
21142115

21152116
// We may have many residuals per shower depending on how many showers ar in the slice
21162117
if (fmShowerRazzle.isValid() && fmShowerRazzle.at(iPart).size()==1) {

‎sbncode/CAFMaker/FillReco.cxx‎

Lines changed: 19 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -4,6 +4,9 @@
44
// \author $Author: [email protected]
55
//////////////////////////////////////////////////////////////////////
66

7+
#include "larcorealg/Geometry/GeometryCore.h"
8+
#include "larcorealg/Geometry/WireReadoutGeom.h"
9+
710
#include "FillReco.h"
811
#include "RecoUtils/RecoUtils.h"
912

@@ -138,6 +141,8 @@ namespace caf
138141
}
139142

140143
void FillCRTPMTMatch(const sbn::crt::CRTPMTMatching &match,
144+
int &topen, int &topex, int &sideen, int &sideex,
145+
//int &bottomen, int &bottomex,
141146
caf::SRCRTPMTMatch &srmatch,
142147
bool allowEmpty){
143148
// allowEmpty does not (yet) matter here
@@ -154,17 +159,24 @@ namespace caf
154159
srmatch.flashPosition = SRVector3D (match.flashPosition.X(), match.flashPosition.Y(), match.flashPosition.Z());
155160
srmatch.flashYWidth = match.flashYWidth;
156161
srmatch.flashZWidth = match.flashZWidth;
157-
srmatch.flashClassification = static_cast<int>(match.flashClassification);
158162
for(const auto& matchedCRTHit : match.matchedCRTHits){
159-
std::cout << "CRTPMTTimeDiff = "<< matchedCRTHit.PMTTimeDiff << "\n";
160163
caf::SRMatchedCRT matchedCRT;
161164
matchedCRT.PMTTimeDiff = matchedCRTHit.PMTTimeDiff;
162165
matchedCRT.time = matchedCRTHit.time;
163166
matchedCRT.sys = matchedCRTHit.sys;
164167
matchedCRT.region = matchedCRTHit.region;
165168
matchedCRT.position = SRVector3D(matchedCRTHit.position.X(), matchedCRTHit.position.Y(), matchedCRTHit.position.Z());
169+
if(matchedCRTHit.PMTTimeDiff < 0){
170+
if(matchedCRTHit.sys == 0) topen++;
171+
else if(matchedCRTHit.sys == 1) sideen++;
172+
}
173+
else if(matchedCRTHit.PMTTimeDiff >= 0){
174+
if(matchedCRTHit.sys == 0) topex++;
175+
else if(matchedCRTHit.sys == 1) sideex++;
176+
}
166177
srmatch.matchedCRTHits.push_back(matchedCRT);
167178
}
179+
srmatch.flashClassification = static_cast<int>(sbn::crt::assignFlashClassification(topen, topex, sideen, sideex, 0, 0));
168180
}
169181

170182

@@ -228,7 +240,7 @@ namespace caf
228240
void FillShowerVars(const recob::Shower& shower,
229241
const recob::Vertex* vertex,
230242
const std::vector<art::Ptr<recob::Hit>> &hits,
231-
const geo::GeometryCore *geom,
243+
const geo::WireReadoutGeom& wireReadout,
232244
unsigned producer,
233245
caf::SRShower &srshower,
234246
bool allowEmpty)
@@ -281,9 +293,9 @@ namespace caf
281293
for(int p = 0; p < 3; ++p) srshower.plane[p].nHits = 0;
282294
for (auto const& hit:hits) ++srshower.plane[hit->WireID().Plane].nHits;
283295

284-
for (geo::PlaneGeo const& plane: geom->Iterate<geo::PlaneGeo>()) {
296+
for (geo::PlaneGeo const& plane: wireReadout.Iterate<geo::PlaneGeo>()) {
285297

286-
const double angleToVert(geom->WireAngleToVertical(plane.View(), plane.ID()) - 0.5*M_PI);
298+
const double angleToVert(wireReadout.WireAngleToVertical(plane.View(), plane.ID()) - 0.5*M_PI);
287299
const double cosgamma(std::abs(std::sin(angleToVert)*shower.Direction().Y()+std::cos(angleToVert)*shower.Direction().Z()));
288300

289301
srshower.plane[plane.ID().Plane].wirePitch = plane.WirePitch()/cosgamma;
@@ -695,7 +707,6 @@ namespace caf
695707
}
696708

697709
void FillTrackChi2PID(const std::vector<art::Ptr<anab::ParticleID>> particleIDs,
698-
const geo::GeometryCore *geom,
699710
caf::SRTrack& srtrack,
700711
bool allowEmpty)
701712
{
@@ -753,7 +764,7 @@ namespace caf
753764
p.wire = h->WireID().Wire;
754765
p.tpc = h->WireID().TPC;
755766
p.channel = h->Channel();
756-
p.sumadc = h->SummedADC();
767+
p.sumadc = h->ROISummedADC();
757768
p.integral = h->Integral();
758769
p.t = h->PeakTime();
759770
p.width = h->RMS();
@@ -815,7 +826,7 @@ namespace caf
815826
void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
816827
const std::vector<art::Ptr<recob::Hit>> &hits,
817828
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
818-
const geo::GeometryCore *geom, const detinfo::DetectorPropertiesData &dprop,
829+
const detinfo::DetectorPropertiesData &dprop,
819830
caf::SRTrack& srtrack,
820831
bool allowEmpty)
821832
{

‎sbncode/CAFMaker/FillReco.h‎

Lines changed: 7 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,3 @@
1-
21
#ifndef CAF_FILLRECO_H
32
#define CAF_FILLRECO_H
43

@@ -9,8 +8,7 @@
98
#include "lardataalg/DetectorInfo/DetectorPropertiesStandard.h"
109

1110
// LArSoft includes
12-
#include "larcore/Geometry/Geometry.h"
13-
#include "larcorealg/Geometry/GeometryCore.h"
11+
#include "larcorealg/Geometry/fwd.h"
1412

1513
#include "lardataobj/RecoBase/PFParticle.h"
1614
#include "lardataobj/RecoBase/Shower.h"
@@ -59,7 +57,7 @@ namespace caf
5957
void FillShowerVars(const recob::Shower& shower,
6058
const recob::Vertex* vertex,
6159
const std::vector<art::Ptr<recob::Hit>> &hits,
62-
const geo::GeometryCore *geom,
60+
const geo::WireReadoutGeom& wireReadout,
6361
unsigned producer,
6462
caf::SRShower& srshower,
6563
bool allowEmpty = false);
@@ -165,7 +163,6 @@ namespace caf
165163

166164
void FillPlaneChi2PID(const anab::ParticleID &particle_id, caf::SRTrkChi2PID &srpid);
167165
void FillTrackChi2PID(const std::vector<art::Ptr<anab::ParticleID>> particleIDs,
168-
const geo::GeometryCore *geom,
169166
caf::SRTrack& srtrack,
170167
bool allowEmpty = false);
171168

@@ -190,7 +187,7 @@ namespace caf
190187
void FillTrackCalo(const std::vector<art::Ptr<anab::Calorimetry>> &calos,
191188
const std::vector<art::Ptr<recob::Hit>> &hits,
192189
bool fill_calo_points, float fillhit_rrstart, float fillhit_rrend,
193-
const geo::GeometryCore *geom, const detinfo::DetectorPropertiesData &dprop,
190+
const detinfo::DetectorPropertiesData &dprop,
194191
caf::SRTrack& srtrack,
195192
bool allowEmpty = false);
196193

@@ -225,8 +222,10 @@ namespace caf
225222
caf::SROpFlash &srflash,
226223
bool allowEmpty = false);
227224
void FillCRTPMTMatch(const sbn::crt::CRTPMTMatching &match,
228-
caf::SRCRTPMTMatch &srmatch,
229-
bool allowEmpty = false);
225+
int &topen, int &topex, int &sideen, int &sidex,
226+
//int &bottomen, int &bottomex,
227+
caf::SRCRTPMTMatch &srmatch,
228+
bool allowEmpty = false);
230229

231230
void FillTPCPMTBarycenterMatch(const sbn::TPCPMTBarycenterMatch *matchInfo,
232231
caf::SRSlice& slice);

‎sbncode/CAFMaker/FillTrue.cxx‎

Lines changed: 15 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1,6 +1,8 @@
11
#include "FillTrue.h"
22

33
#include "larcorealg/GeoAlgo/GeoAlgo.h"
4+
#include "larcorealg/Geometry/GeometryCore.h"
5+
#include "larcorealg/Geometry/WireReadoutGeom.h"
46
#include "larevt/SpaceCharge/SpaceCharge.h"
57
#include "larcore/CoreUtils/ServiceUtil.h"
68
#include "RecoUtils/RecoUtils.h"
@@ -136,7 +138,8 @@ namespace caf {
136138
// Assumes truth matching and calo-points are filled
137139
void FillTrackCaloTruth(const std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> &id_to_ide_map,
138140
const std::vector<simb::MCParticle> &mc_particles,
139-
const geo::GeometryCore *geo,
141+
const geo::GeometryCore& geometry,
142+
const geo::WireReadoutGeom& wireReadout,
140143
const detinfo::DetectorClocksData &clockData,
141144
const spacecharge::SpaceCharge *sce,
142145
caf::SRTrack& srtrack) {
@@ -163,7 +166,7 @@ namespace caf {
163166
const std::vector<std::pair<geo::WireID, const sim::IDE*>> &match_ides = id_to_ide_map.at(srtrack.truth.p.G4ID);
164167
std::map<unsigned, std::vector<const sim::IDE *>> chan_2_ides;
165168
for (auto const &ide_pair: match_ides) {
166-
chan_2_ides[geo->PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second);
169+
chan_2_ides[wireReadout.PlaneWireToChannel(ide_pair.first)].push_back(ide_pair.second);
167170
}
168171

169172
// pre-compute partial ranges
@@ -267,18 +270,19 @@ namespace caf {
267270
// directly apply the true direction. We include the effect of space charge
268271
// on the true pitch here
269272
if (closest_dist >= 0. && direction.Mag() > 1e-4) {
270-
geo::PlaneID plane(srtrack.producer, p.tpc, iplane);
271-
float angletovert = geo->WireAngleToVertical(geo->View(plane), plane) - 0.5*::util::pi<>();
272-
TVector3 loc_mdx_v = loc_nosce_v - direction * (geo->WirePitch(geo->View(plane)) / 2.);
273-
TVector3 loc_pdx_v = loc_nosce_v + direction * (geo->WirePitch(geo->View(plane)) / 2.);
273+
geo::PlaneID planeid(srtrack.producer, p.tpc, iplane);
274+
geo::PlaneGeo const& plane = wireReadout.Plane(planeid);
275+
float angletovert = wireReadout.WireAngleToVertical(plane.View(), planeid) - 0.5*::util::pi<>();
276+
TVector3 loc_mdx_v = loc_nosce_v - direction * (plane.WirePitch() / 2.);
277+
TVector3 loc_pdx_v = loc_nosce_v + direction * (plane.WirePitch() / 2.);
274278

275279
// Convert types for helper functions
276280
geo::Point_t loc_mdx(loc_mdx_v.X(), loc_mdx_v.Y(), loc_mdx_v.Z());
277281
geo::Point_t loc_pdx(loc_pdx_v.X(), loc_pdx_v.Y(), loc_pdx_v.Z());
278282

279283
// Map to wires
280284
if (sce && sce->EnableSimSpatialSCE()) {
281-
int corr = geo->TPC(plane).DriftDir().X();
285+
int corr = geometry.TPC(plane.ID()).DriftDir().X();
282286

283287
geo::Vector_t offset_m = sce->GetPosOffsets(loc_mdx);
284288
offset_m.SetX(offset_m.X()*corr); // convert from drift direction to detector direction
@@ -296,7 +300,7 @@ namespace caf {
296300
double cosgamma = std::abs(std::sin(angletovert)*dir.Y() + std::cos(angletovert)*dir.Z());
297301
double pitch;
298302
if (cosgamma) {
299-
pitch = geo->WirePitch(geo->View(plane))/cosgamma;
303+
pitch = plane.View()/cosgamma;
300304
}
301305
else {
302306
pitch = 0.;
@@ -306,7 +310,7 @@ namespace caf {
306310
geo::Point_t loc_atwires_alongpitch = loc_sce + dir*pitch;
307311
geo::Point_t loc_alongpitch;
308312
if (sce && sce->EnableSimSpatialSCE()) {
309-
geo::Vector_t offset = sce->GetCalPosOffsets(loc_atwires_alongpitch, plane.TPC);
313+
geo::Vector_t offset = sce->GetCalPosOffsets(loc_atwires_alongpitch, planeid.TPC);
310314
loc_alongpitch = loc_atwires_alongpitch + offset;
311315
}
312316
else loc_alongpitch = loc_atwires_alongpitch;
@@ -878,13 +882,13 @@ namespace caf {
878882
return ret;
879883
}
880884

881-
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> PrepSimChannels(const std::vector<art::Ptr<sim::SimChannel>> &simchannels, const geo::GeometryCore &geo) {
885+
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) {
882886
std::map<int, std::vector<std::pair<geo::WireID, const sim::IDE*>>> ret;
883887

884888
for (const art::Ptr<sim::SimChannel> sc : simchannels) {
885889
// Lookup the wire of this channel
886890
raw::ChannelID_t channel = sc->Channel();
887-
std::vector<geo::WireID> maybewire = geo.ChannelToWire(channel);
891+
std::vector<geo::WireID> maybewire = wireReadout.ChannelToWire(channel);
888892
geo::WireID thisWire; // Default constructor makes invalid wire
889893
if (maybewire.size()) thisWire = maybewire[0];
890894

0 commit comments

Comments
 (0)