Skip to content

Commit 5ae81e8

Browse files
nathanielerowenathanielerowe
authored
Merge branch 'develop' into feature/hlay_crt_saturation
2 parents 3e9df79 + f23c4b0 commit 5ae81e8

20 files changed

Lines changed: 953 additions & 469 deletions

‎sbndcode/CRT/CRTAna/CRTAnalysis_module.cc‎

Lines changed: 75 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -33,6 +33,7 @@
3333
#include "sbnobj/SBND/CRT/CRTCluster.hh"
3434
#include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
3535
#include "sbnobj/SBND/CRT/CRTTrack.hh"
36+
#include "sbnobj/SBND/CRT/CRTBlob.hh"
3637
#include "sbnobj/SBND/Timing/DAQTimestamp.hh"
3738

3839
#include "sbndcode/Geometry/GeometryWrappers/CRTGeoService.h"
@@ -87,6 +88,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
8788
void AnalyseCRTTracks(const art::Event &e, const std::vector<art::Ptr<CRTTrack>> &CRTTrackVec, const art::FindManyP<CRTSpacePoint> &tracksToSpacePoints,
8889
const art::FindOneP<CRTCluster> &spacePointsToClusters, const art::FindManyP<CRTStripHit> &clustersToStripHits);
8990

91+
void AnalyseCRTBlobs(std::vector<art::Ptr<CRTBlob>> &CRTBlobVec);
92+
9093
void AnalyseTPCMatching(const art::Event &e, const art::Handle<std::vector<recob::Track>> &TPCTrackHandle,
9194
const art::Handle<std::vector<CRTSpacePoint>> &CRTSpacePointHandle, const art::Handle<std::vector<CRTCluster>> &CRTClusterHandle,
9295
const art::Handle<std::vector<recob::PFParticle>> &PFPHandle,
@@ -104,10 +107,10 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
104107
CRTBackTrackerAlg fCRTBackTrackerAlg;
105108

106109
std::string fMCParticleModuleLabel, fSimDepositModuleLabel, fFEBDataModuleLabel, fCRTStripHitModuleLabel,
107-
fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fTPCTrackModuleLabel,
110+
fCRTClusterModuleLabel, fCRTSpacePointModuleLabel, fCRTTrackModuleLabel, fCRTBlobModuleLabel, fTPCTrackModuleLabel,
108111
fCRTSpacePointMatchingModuleLabel, fCRTTrackMatchingModuleLabel, fPFPModuleLabel, fPTBModuleLabel,
109112
fTDCModuleLabel, fTimingReferenceModuleLabel;
110-
bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fTruthMatch;
113+
bool fDebug, fDataMode, fNoTPC, fHasPTB, fHasTDC, fHasBlobs, fTruthMatch;
111114
//! Adding some of the reco parameters to save corrections
112115
double fPEAttenuation, fTimeWalkNorm, fTimeWalkScale, fPropDelay;
113116

@@ -296,6 +299,16 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
296299
std::vector<double> _tr_truth_theta;
297300
std::vector<double> _tr_truth_phi;
298301

302+
// crt blob information
303+
std::vector<double> _bl_ts0;
304+
std::vector<double> _bl_ets0;
305+
std::vector<double> _bl_ts1;
306+
std::vector<double> _bl_ets1;
307+
std::vector<double> _bl_pe;
308+
std::vector<int> _bl_nsps;
309+
std::vector<std::vector<int>> _bl_nsps_per_tagger;
310+
311+
// tpc track information (including crt matching)
299312
std::vector<double> _tpc_start_x;
300313
std::vector<double> _tpc_start_y;
301314
std::vector<double> _tpc_start_z;
@@ -341,12 +354,14 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
341354
std::vector<double> _tpc_tr_end_z;
342355
std::vector<double> _tpc_tr_score;
343356

357+
// ptb information (trigger board)
344358
std::vector<uint64_t> _ptb_hlt_trigger;
345359
std::vector<uint64_t> _ptb_hlt_timestamp;
346360

347361
std::vector<uint64_t> _ptb_llt_trigger;
348362
std::vector<uint64_t> _ptb_llt_timestamp;
349363

364+
// spec tdc information (timing board)
350365
std::vector<uint32_t> _tdc_channel;
351366
std::vector<uint64_t> _tdc_timestamp;
352367
std::vector<uint64_t> _tdc_offset;
@@ -364,6 +379,7 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
364379
fCRTClusterModuleLabel = p.get<std::string>("CRTClusterModuleLabel", "crtclustering");
365380
fCRTSpacePointModuleLabel = p.get<std::string>("CRTSpacePointModuleLabel", "crtspacepoints");
366381
fCRTTrackModuleLabel = p.get<std::string>("CRTTrackModuleLabel", "crttracks");
382+
fCRTBlobModuleLabel = p.get<std::string>("CRTBlobModuleLabel", "crtblobs");
367383
fTPCTrackModuleLabel = p.get<std::string>("TPCTrackModuleLabel", "pandoraSCETrack");
368384
fCRTSpacePointMatchingModuleLabel = p.get<std::string>("CRTSpacePointMatchingModuleLabel", "crtspacepointmatchingSCE");
369385
fCRTTrackMatchingModuleLabel = p.get<std::string>("CRTTrackMatchingModuleLabel", "crttrackmatchingSCE");
@@ -376,8 +392,8 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
376392
fNoTPC = p.get<bool>("NoTPC", false);
377393
fHasPTB = p.get<bool>("HasPTB", false);
378394
fHasTDC = p.get<bool>("HasTDC", false);
395+
fHasBlobs = p.get<bool>("HasBlobs", false);
379396
fTruthMatch = p.get<bool>("TruthMatch", true);
380-
//! Adding some of the reco parameters to save corrections
381397
fPEAttenuation = p.get<double>("PEAttenuation", 1.0);
382398
fTimeWalkNorm = p.get<double>("TimeWalkNorm", 0.0);
383399
fTimeWalkScale = p.get<double>("TimeWalkScale", 0.0);
@@ -578,6 +594,17 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
578594
fTree->Branch("tr_truth_phi", "std::vector<double>", &_tr_truth_phi);
579595
}
580596

597+
if(fHasBlobs)
598+
{
599+
fTree->Branch("bl_ts0", "std::vector<double>", &_bl_ts0);
600+
fTree->Branch("bl_ets0", "std::vector<double>", &_bl_ets0);
601+
fTree->Branch("bl_ts1", "std::vector<double>", &_bl_ts1);
602+
fTree->Branch("bl_ets1", "std::vector<double>", &_bl_ets1);
603+
fTree->Branch("bl_pe", "std::vector<double>", &_bl_pe);
604+
fTree->Branch("bl_nsps", "std::vector<int>", &_bl_nsps);
605+
fTree->Branch("bl_nsps_per_tagger", "std::vector<std::vector<int>>", &_bl_nsps_per_tagger);
606+
}
607+
581608
if(!fNoTPC)
582609
{
583610
fTree->Branch("tpc_start_x", "std::vector<double>", &_tpc_start_x);
@@ -835,6 +862,23 @@ void sbnd::crt::CRTAnalysis::analyze(art::Event const& e)
835862
// Fill CRTTrack variables
836863
AnalyseCRTTracks(e, CRTTrackVec, tracksToSpacePoints, spacepointsToClusters, clustersToStripHits);
837864

865+
if(fHasBlobs)
866+
{
867+
// Get CRTBlobs
868+
art::Handle<std::vector<CRTBlob>> CRTBlobHandle;
869+
e.getByLabel(fCRTBlobModuleLabel, CRTBlobHandle);
870+
if(!CRTBlobHandle.isValid()){
871+
std::cout << "CRTBlob product " << fCRTBlobModuleLabel << " not found..." << std::endl;
872+
throw std::exception();
873+
}
874+
875+
std::vector<art::Ptr<CRTBlob>> CRTBlobVec;
876+
art::fill_ptr_vector(CRTBlobVec, CRTBlobHandle);
877+
878+
// Fill CRTBlob variables
879+
AnalyseCRTBlobs(CRTBlobVec);
880+
}
881+
838882
if(fNoTPC)
839883
{
840884
// Fill the Tree
@@ -1585,6 +1629,34 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTTracks(const art::Event &e, const std::ve
15851629
}
15861630
}
15871631

1632+
void sbnd::crt::CRTAnalysis::AnalyseCRTBlobs(std::vector<art::Ptr<CRTBlob>> &CRTBlobVec)
1633+
{
1634+
const unsigned nBlobs = CRTBlobVec.size();
1635+
1636+
_bl_ts0.resize(nBlobs);
1637+
_bl_ets0.resize(nBlobs);
1638+
_bl_ts1.resize(nBlobs);
1639+
_bl_ets1.resize(nBlobs);
1640+
_bl_pe.resize(nBlobs);
1641+
_bl_nsps.resize(nBlobs);
1642+
_bl_nsps_per_tagger.resize(nBlobs, std::vector<int>(7));
1643+
1644+
for(unsigned i = 0; i < nBlobs; ++i)
1645+
{
1646+
const auto blob = CRTBlobVec[i];
1647+
1648+
_bl_ts0[i] = blob->Ts0();
1649+
_bl_ets0[i] = blob->Ts0Err();
1650+
_bl_ts1[i] = blob->Ts1();
1651+
_bl_ets1[i] = blob->Ts1Err();
1652+
_bl_pe[i] = blob->PE();
1653+
_bl_nsps[i] = blob->TotalSpacePoints();
1654+
1655+
for(unsigned j = 0; j < 7; ++j)
1656+
_bl_nsps_per_tagger[i][j] = blob->SpacePointsInTagger((CRTTagger)j);
1657+
}
1658+
}
1659+
15881660
void sbnd::crt::CRTAnalysis::AnalyseTPCMatching(const art::Event &e, const art::Handle<std::vector<recob::Track>> &TPCTrackHandle,
15891661
const art::Handle<std::vector<CRTSpacePoint>> &CRTSpacePointHandle, const art::Handle<std::vector<CRTCluster>> &CRTClusterHandle,
15901662
const art::Handle<std::vector<recob::PFParticle>> &PFPHandle,

‎sbndcode/CRT/CRTReco/CMakeLists.txt‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -37,4 +37,10 @@ simple_plugin(
3737
Eigen3::Eigen
3838
)
3939

40+
simple_plugin(
41+
CRTBlobProducer module
42+
lardata::Utilities
43+
sbnobj::SBND_CRT
44+
)
45+
4046
install_fhicl()

‎sbndcode/CRT/CRTReco/CRTBlobProducer_module.cc‎

Lines changed: 219 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,219 @@
1+
////////////////////////////////////////////////////////////////////////
2+
// Class: CRTBlobProducer
3+
// Plugin Type: producer
4+
// File: CRTBlobProducer_module.cc
5+
//
6+
// Author: Henry Lay (h.lay@sheffield.ac.uk)
7+
////////////////////////////////////////////////////////////////////////
8+
9+
#include "art/Framework/Core/EDProducer.h"
10+
#include "art/Framework/Core/ModuleMacros.h"
11+
#include "art/Framework/Principal/Event.h"
12+
#include "art/Framework/Principal/Handle.h"
13+
#include "art/Framework/Principal/Run.h"
14+
#include "art/Framework/Principal/SubRun.h"
15+
#include "canvas/Utilities/InputTag.h"
16+
#include "canvas/Persistency/Common/FindOneP.h"
17+
#include "fhiclcpp/ParameterSet.h"
18+
#include "messagefacility/MessageLogger/MessageLogger.h"
19+
20+
#include "lardata/Utilities/AssociationUtil.h"
21+
22+
#include "sbnobj/SBND/CRT/CRTCluster.hh"
23+
#include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
24+
#include "sbnobj/SBND/CRT/CRTBlob.hh"
25+
26+
#include <numeric>
27+
28+
namespace sbnd::crt {
29+
class CRTBlobProducer;
30+
}
31+
32+
class sbnd::crt::CRTBlobProducer : public art::EDProducer {
33+
public:
34+
explicit CRTBlobProducer(fhicl::ParameterSet const& p);
35+
36+
CRTBlobProducer(CRTBlobProducer const&) = delete;
37+
CRTBlobProducer(CRTBlobProducer&&) = delete;
38+
CRTBlobProducer& operator=(CRTBlobProducer const&) = delete;
39+
CRTBlobProducer& operator=(CRTBlobProducer&&) = delete;
40+
41+
void produce(art::Event& e) override;
42+
43+
void OrderSpacePoints(std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec, const art::FindOneP<CRTCluster> &spacePointsToCluster);
44+
45+
std::vector<std::pair<CRTBlob, std::set<unsigned>>> CreateBlobCandidates(const std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec,
46+
const art::FindOneP<CRTCluster> &spacePointsToCluster);
47+
48+
void TimeErrorCalculator(const std::vector<double> &times, double &mean, double &err);
49+
50+
void OrderBlobCandidates(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates);
51+
52+
std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> ChoseBlobs(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates);
53+
54+
private:
55+
56+
std::string fCRTSpacePointModuleLabel;
57+
double fCoincidenceTimeRequirement;
58+
bool fUseTs0;
59+
};
60+
61+
62+
sbnd::crt::CRTBlobProducer::CRTBlobProducer(fhicl::ParameterSet const& p)
63+
: EDProducer{p}
64+
, fCRTSpacePointModuleLabel(p.get<std::string>("CRTSpacePointModuleLabel"))
65+
, fCoincidenceTimeRequirement(p.get<double>("CoincidenceTimeRequirement"))
66+
, fUseTs0(p.get<bool>("UseTs0"))
67+
{
68+
produces<std::vector<CRTBlob>>();
69+
produces<art::Assns<CRTSpacePoint, CRTBlob>>();
70+
}
71+
72+
void sbnd::crt::CRTBlobProducer::produce(art::Event& e)
73+
{
74+
auto blobVec = std::make_unique<std::vector<CRTBlob>>();
75+
auto blobSpacePointAssn = std::make_unique<art::Assns<CRTSpacePoint, CRTBlob>>();
76+
77+
art::Handle<std::vector<CRTSpacePoint>> CRTSpacePointHandle;
78+
e.getByLabel(fCRTSpacePointModuleLabel, CRTSpacePointHandle);
79+
80+
std::vector<art::Ptr<CRTSpacePoint>> CRTSpacePointVec;
81+
art::fill_ptr_vector(CRTSpacePointVec, CRTSpacePointHandle);
82+
83+
art::FindOneP<CRTCluster> spacePointsToCluster(CRTSpacePointHandle, e, fCRTSpacePointModuleLabel);
84+
85+
OrderSpacePoints(CRTSpacePointVec, spacePointsToCluster);
86+
87+
std::vector<std::pair<CRTBlob, std::set<unsigned>>> blobCandidates = CreateBlobCandidates(CRTSpacePointVec, spacePointsToCluster);
88+
89+
OrderBlobCandidates(blobCandidates);
90+
91+
std::vector<std::pair<CRTBlob, std::set<unsigned>>> chosenBlobs = ChoseBlobs(blobCandidates);
92+
93+
for(auto const& [blob, spIDs] : chosenBlobs)
94+
{
95+
blobVec->push_back(blob);
96+
97+
for(auto const& spID : spIDs)
98+
util::CreateAssn(*this, e, *blobVec, CRTSpacePointVec[spID], *blobSpacePointAssn);
99+
}
100+
101+
e.put(std::move(blobVec));
102+
e.put(std::move(blobSpacePointAssn));
103+
}
104+
105+
void sbnd::crt::CRTBlobProducer::OrderSpacePoints(std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec, const art::FindOneP<CRTCluster> &spacePointsToCluster)
106+
{
107+
std::sort(spacePointVec.begin(), spacePointVec.end(),
108+
[&](const art::Ptr<CRTSpacePoint> &a, const art::Ptr<CRTSpacePoint> &b) -> bool {
109+
if(fUseTs0)
110+
return a->Ts0() < b->Ts0();
111+
else
112+
return a->Ts1() < b->Ts1();
113+
});
114+
}
115+
116+
std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> sbnd::crt::CRTBlobProducer::CreateBlobCandidates(const std::vector<art::Ptr<CRTSpacePoint>> &spacePointVec,
117+
const art::FindOneP<CRTCluster> &spacePointsToCluster)
118+
{
119+
std::vector<std::pair<CRTBlob, std::set<unsigned>>> candidates;
120+
121+
for(unsigned i = 0; i < spacePointVec.size(); ++i)
122+
{
123+
const art::Ptr<CRTSpacePoint> primarySpacePoint = spacePointVec[i];
124+
const art::Ptr<CRTCluster> primaryCluster = spacePointsToCluster.at(primarySpacePoint.key());
125+
126+
std::set<unsigned> used_spacepoints = { i };
127+
std::vector<double> t0s = { primarySpacePoint->Ts0() };
128+
std::vector<double> t1s = { primarySpacePoint->Ts1() };
129+
std::vector<double> pes = { primarySpacePoint->PE() };
130+
131+
std::map<CRTTagger, int> taggers;
132+
for(int i = 0; i < 7; ++i)
133+
taggers[(CRTTagger)i] = 0;
134+
++taggers[primaryCluster->Tagger()];
135+
136+
for(unsigned ii = i+1; ii < spacePointVec.size(); ++ii)
137+
{
138+
const art::Ptr<CRTSpacePoint> secondarySpacePoint = spacePointVec[ii];
139+
const art::Ptr<CRTCluster> secondaryCluster = spacePointsToCluster.at(secondarySpacePoint.key());
140+
141+
const double tdiff_prim_sec = fUseTs0 ? secondarySpacePoint->Ts0() - primarySpacePoint->Ts0() : secondarySpacePoint->Ts1() - primarySpacePoint->Ts1();
142+
143+
if(tdiff_prim_sec > fCoincidenceTimeRequirement)
144+
break;
145+
146+
used_spacepoints.insert(ii);
147+
t0s.push_back(secondarySpacePoint->Ts0());
148+
t1s.push_back(secondarySpacePoint->Ts1());
149+
pes.push_back(secondarySpacePoint->PE());
150+
++taggers[secondaryCluster->Tagger()];
151+
}
152+
153+
double t0, et0;
154+
TimeErrorCalculator(t0s, t0, et0);
155+
156+
double t1, et1;
157+
TimeErrorCalculator(t1s, t1, et1);
158+
159+
const double pe = std::accumulate(pes.begin(), pes.end(), 0.);
160+
161+
const CRTBlob blob(t0, et0, t1, et1, pe, taggers);
162+
163+
candidates.emplace_back(blob, used_spacepoints);
164+
}
165+
166+
return candidates;
167+
}
168+
169+
void sbnd::crt::CRTBlobProducer::TimeErrorCalculator(const std::vector<double> &times, double &mean, double &err)
170+
{
171+
double sum = 0.;
172+
for(auto const &time : times)
173+
sum += time;
174+
175+
mean = sum / times.size();
176+
177+
double summed_var = 0.;
178+
for(auto const &time : times)
179+
summed_var += std::pow((time - mean), 2);
180+
181+
err = std::sqrt(summed_var / (times.size() - 1));
182+
}
183+
184+
void sbnd::crt::CRTBlobProducer::OrderBlobCandidates(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates)
185+
{
186+
std::sort(blobCandidates.begin(), blobCandidates.end(),
187+
[&](const std::pair<CRTBlob, std::set<unsigned>> &a, const std::pair<CRTBlob, std::set<unsigned>> &b) -> bool {
188+
return fUseTs0 ? a.first.Ts0Err() < b.first.Ts0Err() : a.first.Ts1Err() < b.first.Ts1Err();
189+
});
190+
}
191+
192+
std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> sbnd::crt::CRTBlobProducer::ChoseBlobs(std::vector<std::pair<CRTBlob, std::set<unsigned>>> &blobCandidates)
193+
{
194+
std::vector<std::pair<sbnd::crt::CRTBlob, std::set<unsigned>>> chosenBlobs;
195+
196+
std::set<unsigned> used;
197+
198+
for(auto const& [blob, spIDs] : blobCandidates)
199+
{
200+
bool keep = true;
201+
for(auto const& spID : spIDs)
202+
{
203+
if(used.count(spID) != 0)
204+
keep = false;
205+
}
206+
207+
if(keep)
208+
{
209+
chosenBlobs.emplace_back(blob, spIDs);
210+
211+
for(auto const& spID : spIDs)
212+
used.insert(spID);
213+
}
214+
}
215+
216+
return chosenBlobs;
217+
}
218+
219+
DEFINE_ART_MODULE(sbnd::crt::CRTBlobProducer)

0 commit comments

Comments
 (0)