Skip to content

Commit 1eabcf6

Browse files
nathanielerowenathanielerowe
authored
Merge pull request #901 from SBNSoftware/bugfix/BlipsOnEmptyEvents
Bugfix/blips on empty events
2 parents d782c24 + dc7d054 commit 1eabcf6

6 files changed

Lines changed: 135 additions & 142 deletions

File tree

‎sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc‎

Lines changed: 64 additions & 89 deletions
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,6 @@ namespace blip {
3232

3333

3434
// Loop over cryostats
35-
std::cout<<"NCryostats: "<<fGeom.Ncryostats()<<"\n";
3635
for(size_t cstat=0; cstat<fGeom.Ncryostats(); cstat++){
3736
auto const& cryoid = geo::CryostatID(cstat);
3837

@@ -41,16 +40,14 @@ namespace blip {
4140
auto const& tpcid = geo::TPCID(cryoid,tpc);
4241

4342
// Loop planes in TPC 'tpc'
44-
auto const& plane0id = geo::PlaneID(cstat,tpc,0);
45-
auto const& plane0geo = wireReadoutGeom->Get().Plane(plane0id);
43+
auto const& plane0id = geo::PlaneID(cstat,tpc,0);
44+
auto const& plane0geo = wireReadoutGeom->Get().Plane(plane0id);
4645
for(size_t pl=0; pl<wireReadoutGeom->Get().Nplanes(tpcid); pl++){
4746
auto const& planeid = geo::PlaneID(cstat,tpc,pl);
4847
auto const& planegeo = wireReadoutGeom->Get().Plane(planeid);
4948
kNumChannels += planegeo.Nwires();
5049

5150
float offset = detProp.GetXTicksOffset(pl,tpc,cstat);
52-
std::cout<<"CRYOSTAT "<<cstat<<" / TPC "<<tpc<<" / PLANE "<<pl<<": "<<planegeo.Nwires()<<" wires\n";
53-
std::cout<<" XTicksOffset (from detProp): "<<offset<<"\n";
5451

5552
kXTicksOffsets[cstat][tpc][pl] = 0;
5653

@@ -66,12 +63,9 @@ namespace blip {
6663
auto const& tpcgeom = cryostat.TPC(tpc);
6764
auto const xyz = plane0geo.GetCenter();
6865
const double dir((tpcgeom.DriftSign() == geo::DriftSign::Negative) ? +1.0 : -1.0);
69-
float x_ticks_coefficient = kDriftVelocity*kTickPeriod;
66+
float x_ticks_coefficient = kDriftVelocity*kTickPeriod;
7067
float goofy_offset = -xyz.X() / (dir * x_ticks_coefficient);
71-
std::cout<<" After geometric correction: "<<offset - goofy_offset<<"\n";
72-
7368
kXTicksOffsets[cstat][tpc][pl] = offset - goofy_offset;
74-
7569
} else {
7670

7771
// for the case of 2D wirecell workflow, the plane-to-plane
@@ -82,12 +76,7 @@ namespace blip {
8276
// we still want to correct for global trigger offset though:
8377
kXTicksOffsets[cstat][tpc][pl] = clockData.TriggerTime()/kTickPeriod;
8478

85-
}
86-
87-
// additional ad-hoc corrections supplied by user
88-
//kXTicksOffsets[cstat][tpc][pl] += fTimeOffset[pl];
89-
std::cout << " offsetting plane " << pl << " by " << fTimeOffset[pl] << " ticks " << std::endl;
90-
79+
}
9180
}
9281
}
9382
}
@@ -164,18 +153,18 @@ namespace blip {
164153
for(int i=0; i<kNplanes; i++) {
165154
if( i == fCaloPlane ) continue;
166155
h_clust_overlap[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_overlap",iTPC,i), Form("TPC %i, Plane %i clusters;Overlap fraction",iTPC,i),101,0,1.01);
167-
h_clust_dt[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_dt",iTPC,i), Form("TPC %i, Plane %i clusters;dT [ticks]",iTPC,i),200,-10,10);
156+
h_clust_dt[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_dt",iTPC,i), Form("TPC %i, Plane %i clusters;dT [ticks]",iTPC,i),200,-10,10);
168157
h_clust_dtfrac[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_dtfrac",iTPC,i), Form("TPC %i, Plane %i clusters;Charge-weighted mean dT/RMS",iTPC,i),150,-1.5,1.5);
169158

170159
h_clust_q[iTPC][i] = hdir.make<TH2D>(Form("t%i_p%i_clust_charge",iTPC,i),
171160
Form("Pre-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3} e-]",iTPC,fCaloPlane,i),
172161
qbins,0,qmax,qbins,0,qmax);
173-
h_clust_q[iTPC][i]->SetOption("colz");
162+
h_clust_q[iTPC][i]->SetOption("colz");
174163

175164
h_clust_q_cut[iTPC][i] = hdir.make<TH2D>(Form("t%i_p%i_clust_charge_cut",iTPC,i),
176165
Form("Post-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3}]",iTPC,fCaloPlane,i),
177166
qbins,0,qmax,qbins,0,qmax);
178-
h_clust_q_cut[iTPC][i]->SetOption("colz");
167+
h_clust_q_cut[iTPC][i]->SetOption("colz");
179168

180169
h_clust_score[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_matchscore",iTPC,i), Form("TPC %i, Plane %i clusters;Match score",iTPC,i),101,0,1.01);
181170

@@ -189,7 +178,7 @@ namespace blip {
189178
h_clust_truematch_q[iTPC][i] = hdir.make<TH2D>(Form("t%i_p%i_clust_truematch_charge",iTPC,i),
190179
Form("Pre-cut, TPC %i;Plane %i cluster charge [#times10^{3} e-];Plane %i cluster charge [#times10^{3} e-]",iTPC,fCaloPlane,i),
191180
qbins,0,qmax,qbins,0,qmax);
192-
h_clust_truematch_q[iTPC][i]->SetOption("colz");
181+
h_clust_truematch_q[iTPC][i]->SetOption("colz");
193182

194183
h_clust_truematch_score[iTPC][i] = hdir.make<TH1D>(Form("t%i_p%i_clust_truematch_matchscore",iTPC,i), Form("TPC %i, Plane %i clusters;Match score",iTPC,i),101,0,1.01);
195184

@@ -239,6 +228,7 @@ namespace blip {
239228
void BlipRecoAlg::reconfigure( fhicl::ParameterSet const& pset ){
240229

241230
fHitProducer = pset.get<std::string> ("HitProducer", "gaushit");
231+
fHitTruthMatcher = pset.get<std::string> ("HitTruthMatcher", "blipgaushitTruthMatch");
242232
fTrkProducer = pset.get<std::string> ("TrkProducer", "pandora");
243233
fGeantProducer = pset.get<std::string> ("GeantProducer", "largeant");
244234
fSimDepProducer = pset.get<std::string> ("SimEDepProducer", "ionandscint");
@@ -357,45 +347,39 @@ namespace blip {
357347
if (evt.getByLabel(fGeantProducer,pHandle))
358348
art::fill_ptr_vector(plist, pHandle);
359349

360-
// -- SimEnergyDeposits (usually dropped in reco
361-
//art::Handle<std::vector<sim::EnergyDeposit> > sedHandle;
362-
std::vector<sim::IDE > sedlist;
363-
//if (evt.getByLabel(fSimDepProducer,sedHandle)){
364-
// art::fill_ptr_vector(sedlist, sedHandle);
365-
// }
366-
// -- SimChannels (usually dropped in reco)
350+
// -- SimEnergyDeposits
351+
art::Handle<std::vector<sim::SimEnergyDeposit> > sedHandle;
352+
std::vector<art::Ptr<sim::SimEnergyDeposit> > sedlist;
353+
if (evt.getByLabel(fSimDepProducer,sedHandle)){
354+
art::fill_ptr_vector(sedlist, sedHandle);
355+
}
356+
std::vector<sim::IDE > sIDElist;
357+
// -- SimChannels
367358
art::Handle<std::vector<sim::SimChannel> > simchanHandle;
368359
std::vector<art::Ptr<sim::SimChannel> > simchanlist;
369360
if (evt.getByLabel(fSimChanProducer,simchanHandle))
370-
{
361+
{
371362
art::fill_ptr_vector(simchanlist, simchanHandle);
372-
//Loop over channels to get full sedlist
363+
//Loop over channels to get full sIDElist
373364
for(int chIndex=0; chIndex<int(simchanlist.size()); chIndex++)
374-
{
375-
std::vector<geo::WireID> wids = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
376-
const geo::PlaneID& planeID = wids[0].planeID();
377-
if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values
378-
std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, 999999999);
379-
for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
380365
{
381-
//art::fill_ptr_vector(sedlist, simchanHandle.TrackIDsAndEnergies(0, 99999999));
382-
sedlist.push_back( TempChIDE[ideIndex] ); //may need to add a &
366+
std::vector<geo::WireID> wids = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
367+
const geo::PlaneID& planeID = wids[0].planeID();
368+
if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values
369+
unsigned int MaxTDCTick = UINT_MAX;
370+
std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, MaxTDCTick);
371+
for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
372+
{
373+
sIDElist.push_back( TempChIDE[ideIndex] ); //may need to add a &
374+
}
383375
}
384-
}
385-
}
386-
376+
}
387377
// -- hits (from input module, usually track-masked subset of gaushit)
388378
art::Handle< std::vector<recob::Hit> > hitHandle;
389379
std::vector<art::Ptr<recob::Hit> > hitlist;
390380
if (evt.getByLabel(fHitProducer,hitHandle))
391381
art::fill_ptr_vector(hitlist, hitHandle);
392382

393-
// -- hits (from gaushit), these are used in truth-matching of hits
394-
art::Handle< std::vector<recob::Hit> > hitHandleGH;
395-
std::vector<art::Ptr<recob::Hit> > hitlistGH;
396-
if (evt.getByLabel("gaushit",hitHandleGH))
397-
art::fill_ptr_vector(hitlistGH, hitHandleGH);
398-
399383
// -- tracks
400384
art::Handle< std::vector<recob::Track> > tracklistHandle;
401385
std::vector<art::Ptr<recob::Track> > tracklist;
@@ -404,8 +388,8 @@ namespace blip {
404388

405389
// -- associations
406390
art::FindManyP<recob::Track> fmtrk(hitHandle,evt,fTrkProducer);
407-
art::FindManyP<recob::Track> fmtrkGH(hitHandleGH,evt,fTrkProducer);
408-
art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleGH,evt,"gaushitTruthMatch");
391+
art::FindManyP<recob::Track> fmtrkGH(hitHandle,evt,fTrkProducer);
392+
art::FindMany<simb::MCParticle, anab::BackTrackerHitMatchingData> fmhh(hitHandle,evt,fHitTruthMatcher);
409393
/*
410394
//====================================================
411395
// Update map of bad channels for this event
@@ -441,23 +425,7 @@ namespace blip {
441425
// use gaushitTruthMatch later on)
442426
//===============================================================
443427
std::map< int, int > map_gh;
444-
// if input collection is already gaushit, this is trivial
445-
if( fHitProducer == "gaushit" ) {
446-
for(auto& h : hitlist ) map_gh[h.key()] = h.key();
447-
// ... but if not, find the matching gaushit. There's no convenient
448-
// hit ID, so we must loop through and compare channel/time (ugh)
449-
} else {
450-
std::map<int,std::vector<int>> map_chan_ghid;
451-
for(auto& gh : hitlistGH ) map_chan_ghid[gh->Channel()].push_back(gh.key());
452-
for(auto& h : hitlist ) {
453-
for(auto& igh : map_chan_ghid[h->Channel()]){
454-
if( hitlistGH[igh]->PeakTime() != h->PeakTime() ) continue;
455-
map_gh[h.key()] = igh;
456-
break;
457-
}
458-
}
459-
}
460-
428+
for(auto& h : hitlist ) map_gh[h.key()] = h.key();
461429
//=====================================================
462430
// Record PDG for every G4 Track ID
463431
//=====================================================
@@ -526,15 +494,23 @@ namespace blip {
526494
if( plist.size() ) {
527495
pinfo.resize(plist.size());
528496
for(size_t i = 0; i<plist.size(); i++){
529-
BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
497+
//use sim::EnergyDeposits by default. This is heavy and may be dropped
498+
if(sedlist.size()>0)
499+
{
500+
BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
501+
}
502+
else //use sim::Channel -> IDE otherwise. This is usually kept but results in strange bugs.
503+
{
504+
BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sIDElist, fCaloPlane);
505+
}
530506
if( map_g4trkid_charge[pinfo[i].trackId] ) pinfo[i].numElectrons = (int)map_g4trkid_charge[pinfo[i].trackId];
507+
else pinfo[i].numElectrons = pinfo[i].depElectrons; // JACOB ADDITION Jan22, 2026 for strange channel IDE results
531508
pinfo[i].index = i;
532509
}
533510
BlipUtils::MakeTrueBlips(pinfo, trueblips);
534511
BlipUtils::MergeTrueBlips(trueblips, fTrueBlipMergeDist);
535512
}
536513

537-
538514
//=======================================
539515
// Map track IDs to the index in the vector
540516
//=======================================
@@ -647,14 +623,8 @@ namespace blip {
647623

648624

649625
// find associated track
650-
if( fHitProducer == "gaushit" && fmtrk.isValid() ) {
626+
if( fmtrk.isValid() ) {
651627
if(fmtrk.at(i).size()) hitinfo[i].trkid = fmtrk.at(i)[0]->ID();
652-
653-
// if the hit collection didn't have associations made
654-
// to the tracks, try gaushit instead
655-
} else if ( fmtrkGH.isValid() && map_gh.size() ) {
656-
int gi = map_gh[i];
657-
if (fmtrkGH.at(gi).size()) hitinfo[i].trkid= fmtrkGH.at(gi)[0]->ID();
658628
}
659629

660630
// add to the map
@@ -664,8 +634,7 @@ namespace blip {
664634
if( hitinfo[i].trkid < 0 ) nhits_untracked++;
665635
//printf(" %lu plane: %i, wire: %i, time: %i\n",i,hitinfo[i].plane,hitinfo[i].wire,int(hitinfo[i].driftTime));
666636

667-
}//endloop over hits
668-
637+
}//endloop over hits
669638
//for(auto& a : tpc_plane_hitsMap ) {
670639
//for(auto& b : a.second )
671640
//std::cout<<"TPC "<<a.first<<", plane "<<b.first<<": "<<b.second.size()<<" hits\n";
@@ -690,6 +659,7 @@ namespace blip {
690659

691660

692661
// Basic track inclusion cut: exclude hits that were tracked
662+
int Tracked=0;
693663
for(size_t i=0; i<hitlist.size(); i++){
694664
if( hitinfo[i].trkid < 0 ) continue;
695665
auto it = map_trkid_index.find(hitinfo[i].trkid);
@@ -698,9 +668,9 @@ namespace blip {
698668
if( tracklist[trkindex]->Length() > fMaxHitTrkLength ) {
699669
hitIsTracked[i] = true;
700670
hitIsGood[i] = false;
671+
Tracked++;
701672
}
702-
}
703-
673+
}
704674

705675
// Filter based on hit properties. For hits that are a part of
706676
// multi-gaussian fits (multiplicity > 1), need to re-think this.
@@ -731,7 +701,6 @@ namespace blip {
731701
// Hit clustering
732702
// ---------------------------------------------------
733703
std::map<int,std::map<int,std::vector<int>>> tpc_planeclustsMap;
734-
735704
for(auto const& tpc_plane_hitsMap : cryo_tpc_plane_hitsMap ) {
736705

737706
for(auto const& plane_hitsMap : tpc_plane_hitsMap.second ) {
@@ -892,7 +861,6 @@ namespace blip {
892861

893862
float _matchQDiffLimit= (fMatchQDiffLimit <= 0 ) ? std::numeric_limits<float>::max() : fMatchQDiffLimit;
894863
float _matchMaxQRatio = (fMatchMaxQRatio <= 0 ) ? std::numeric_limits<float>::max() : fMatchMaxQRatio;
895-
896864
for(auto& tpcMap : tpc_planeclustsMap ) { // loop on TPCs
897865

898866
//std::cout
@@ -935,15 +903,15 @@ namespace blip {
935903
// Check that the two central wires intersect
936904
// *******************************************
937905
double y, z;
938-
geo::Point_t intsec_p;
939-
std::vector<geo::WireID> A_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcA.CenterChan);
940-
std::vector<geo::WireID> B_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcB.CenterChan);
906+
geo::Point_t intsec_p;
907+
std::vector<geo::WireID> A_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcA.CenterChan);
908+
std::vector<geo::WireID> B_wireids = wireReadoutGeom->Get().ChannelToWire((unsigned int)hcB.CenterChan);
941909

942-
if( !wireReadoutGeom->Get().WireIDsIntersect(A_wireids.at(0),B_wireids.at(0),intsec_p)) continue;
943-
// Save intersect location, so we don't have to
910+
if( !wireReadoutGeom->Get().WireIDsIntersect(A_wireids.at(0),B_wireids.at(0),intsec_p)) continue;
911+
// Save intersect location, so we don't have to
944912
// make another call to the Geometry service later
945-
y = intsec_p.Y();
946-
z = intsec_p.Z();
913+
y = intsec_p.Y();
914+
z = intsec_p.Z();
947915
TVector3 xloc(0,y,z);
948916
hcA.IntersectLocations[hcB.ID] = xloc;
949917
hcB.IntersectLocations[hcA.ID] = xloc;
@@ -1046,7 +1014,6 @@ namespace blip {
10461014
for(auto& hc : hcGroup ) {
10471015
hitclust[hc.ID].isMatched = true;
10481016
for(auto hit : hitclust[hc.ID].HitIDs) hitinfo[hit].ismatch = true;
1049-
10501017
// Diagnostic plots for successful 3-plane matches
10511018
//if( picky && hc.Plane != fCaloPlane ) {
10521019
//float q1 = (float)newBlip.clusters[fCaloPlane].Charge;
@@ -1096,18 +1063,26 @@ namespace blip {
10961063
// if we made it this far, the blip is good!
10971064
// associate this blip with the hits and clusters within it
10981065
newBlip.ID = blips.size();
1099-
blips.push_back(newBlip);
11001066
for(auto& hc : hcGroup ) {
11011067
hitclust[hc.ID].BlipID = newBlip.ID;
11021068
for( auto& h : hc.HitIDs ) hitinfo[h].blipid = newBlip.ID;
11031069
}
1070+
//BLIPS HAVE A COPY OF HITCLUSTERS NOT A POINTER
1071+
//UPDATE THE HITCLUSTER VARS THAT HAVE CHANGED SINCE CONSTRUCTION
1072+
for(int iclust=0; iclust<int(sizeof(newBlip.clusters)/sizeof(newBlip.clusters[0])); iclust++)
1073+
{
1074+
if(newBlip.clusters[iclust].ID>-1){
1075+
newBlip.clusters[iclust].BlipID = newBlip.ID;
1076+
newBlip.clusters[iclust].isMatched = true;
1077+
}
1078+
}
1079+
blips.push_back(newBlip);
11041080

11051081

11061082
}//endif ncands > 0
11071083
}//endloop over caloplane ("Plane A") clusters
11081084
}//endif calo plane has clusters
11091085
}//endloop over TPCs
1110-
11111086
// Re-index the clusters after removing unmatched
11121087
if( !keepAllClusts ) {
11131088
std::vector<blip::HitClust> hitclust_filt;

‎sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -129,6 +129,7 @@ namespace blip {
129129

130130
// --- FCL configs ---
131131
std::string fHitProducer;
132+
std::string fHitTruthMatcher;
132133
std::string fTrkProducer;
133134
std::string fGeantProducer;
134135
std::string fSimDepProducer;

0 commit comments

Comments
 (0)