Bugfix/blips on empty events

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc

@@ -32,7 +32,6 @@ namespace blip {
 
 
     // Loop over cryostats
-    std::cout<<"NCryostats: "<<fGeom.Ncryostats()<<"\n";
     for(size_t cstat=0; cstat<fGeom.Ncryostats(); cstat++){
       auto const& cryoid = geo::CryostatID(cstat);
 
@@ -49,8 +48,6 @@ namespace blip {
           kNumChannels += planegeo.Nwires();
 
           float offset = detProp.GetXTicksOffset(pl,tpc,cstat);
-          std::cout<<"CRYOSTAT "<<cstat<<" / TPC "<<tpc<<" / PLANE "<<pl<<":  "<<planegeo.Nwires()<<" wires\n";
-          std::cout<<"  XTicksOffset (from detProp): "<<offset<<"\n";
 
           kXTicksOffsets[cstat][tpc][pl] = 0;
 
@@ -68,10 +65,7 @@ namespace blip {
             const double dir((tpcgeom.DriftSign() == geo::DriftSign::Negative) ? +1.0 : -1.0);
 	    float x_ticks_coefficient = kDriftVelocity*kTickPeriod;
             float goofy_offset = -xyz.X() / (dir * x_ticks_coefficient);
-	    std::cout<<"  After geometric correction: "<<offset - goofy_offset<<"\n";
-
             kXTicksOffsets[cstat][tpc][pl] = offset - goofy_offset;
-
 	  } else {
 
             // for the case of 2D wirecell workflow, the plane-to-plane
@@ -85,9 +79,7 @@ namespace blip {
           }
 
           // additional ad-hoc corrections supplied by user
-          //kXTicksOffsets[cstat][tpc][pl] += fTimeOffset[pl];
-          std::cout << " offsetting plane " << pl << " by " << fTimeOffset[pl] << " ticks " << std::endl;
-
+          //kXTicksOffsets[cstat][tpc][pl] += fTimeOffset[pl];          
         }
       }
     }
@@ -357,33 +349,33 @@ namespace blip {
     if (evt.getByLabel(fGeantProducer,pHandle))
       art::fill_ptr_vector(plist, pHandle);
 
-    // -- SimEnergyDeposits (usually dropped in reco
-    //art::Handle<std::vector<sim::EnergyDeposit> > sedHandle;
-    std::vector<sim::IDE > sedlist;
-    //if (evt.getByLabel(fSimDepProducer,sedHandle)){
-    //  art::fill_ptr_vector(sedlist, sedHandle);
-    // }
-    // -- SimChannels (usually dropped in reco)
+    // -- SimEnergyDeposits 
+    art::Handle<std::vector<sim::SimEnergyDeposit> > sedHandle;
+    std::vector<art::Ptr<sim::SimEnergyDeposit> > sedlist;
+    if (evt.getByLabel(fSimDepProducer,sedHandle)){
+      art::fill_ptr_vector(sedlist, sedHandle);
+     }
+    std::vector<sim::IDE > sIDElist;
+    // -- SimChannels 
     art::Handle<std::vector<sim::SimChannel> > simchanHandle;
     std::vector<art::Ptr<sim::SimChannel> > simchanlist;
     if (evt.getByLabel(fSimChanProducer,simchanHandle))
-      { 
+    { 
       art::fill_ptr_vector(simchanlist, simchanHandle);
-      //Loop over channels to get full sedlist
+      //Loop over channels to get full sIDElist
       for(int chIndex=0; chIndex<int(simchanlist.size()); chIndex++)
-	{
-	  std::vector<geo::WireID> wids    = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
-	  const geo::PlaneID&      planeID = wids[0].planeID();
-	  if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values 
-	  std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, 999999999);
-	  for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
 	    {
-	      //art::fill_ptr_vector(sedlist, simchanHandle.TrackIDsAndEnergies(0, 99999999));
-	      sedlist.push_back( TempChIDE[ideIndex] ); //may need to add a &
+      std::vector<geo::WireID> wids    = wireReadoutGeom->Get().ChannelToWire( (*(simchanlist[chIndex])).Channel() ); //Not sure why this is a vector, but it should have len 1
+      const geo::PlaneID&      planeID = wids[0].planeID();
+      if(int(planeID.Plane) != fCaloPlane) continue; //only take calorimetry plane IDE values 
+      unsigned int MaxTDCTick = 4294967294; // 1 under 32 bit unsigned int max
+      std::vector< sim::IDE > TempChIDE = (*simchanlist[chIndex]).TrackIDsAndEnergies(0, MaxTDCTick);
+      for(int ideIndex=0; ideIndex<int(TempChIDE.size()); ideIndex++)
+	        {
+	        sIDElist.push_back( TempChIDE[ideIndex] ); //may need to add a &
+	       }
 	    }
-	}
-      }
-
+    }
     // -- hits (from input module, usually track-masked subset of gaushit)
     art::Handle< std::vector<recob::Hit> > hitHandle;
     std::vector<art::Ptr<recob::Hit> > hitlist;
@@ -393,8 +385,8 @@ namespace blip {
     // -- hits (from gaushit), these are used in truth-matching of hits
     art::Handle< std::vector<recob::Hit> > hitHandleGH;
     std::vector<art::Ptr<recob::Hit> > hitlistGH;
-    if (evt.getByLabel("gaushit",hitHandleGH))
-      art::fill_ptr_vector(hitlistGH, hitHandleGH);
+    if(evt.getByLabel("specialblipgaushit",hitHandleGH)) art::fill_ptr_vector(hitlistGH, hitHandleGH);
+    else if(evt.getByLabel("gaushit",hitHandleGH))       art::fill_ptr_vector(hitlistGH, hitHandleGH);
 
     // -- tracks
     art::Handle< std::vector<recob::Track> > tracklistHandle;
@@ -405,7 +397,7 @@ namespace blip {
     // -- associations
     art::FindManyP<recob::Track> fmtrk(hitHandle,evt,fTrkProducer);
     art::FindManyP<recob::Track> fmtrkGH(hitHandleGH,evt,fTrkProducer);
-    art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleGH,evt,"gaushitTruthMatch");
+    art::FindMany<simb::MCParticle,anab::BackTrackerHitMatchingData> fmhh(hitHandleGH,evt,"blipgaushitTruthMatch");
     /*
     //====================================================
     // Update map of bad channels for this event
@@ -442,7 +434,7 @@ namespace blip {
     //===============================================================
     std::map< int, int > map_gh;
     // if input collection is already gaushit, this is trivial
-    if( fHitProducer == "gaushit" ) {
+    if( fHitProducer == "gaushit" ||  fHitProducer == "specialblipgaushit") {
       for(auto& h : hitlist ) map_gh[h.key()] = h.key(); 
     // ... but if not, find the matching gaushit. There's no convenient
     // hit ID, so we must loop through and compare channel/time (ugh)
@@ -456,8 +448,7 @@ namespace blip {
           break;
         }
       }
-    }
-
+    }   
     //=====================================================
     // Record PDG for every G4 Track ID
     //=====================================================
@@ -526,15 +517,25 @@ namespace blip {
     if( plist.size() ) {
       pinfo.resize(plist.size());
       for(size_t i = 0; i<plist.size(); i++){
-        BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
+        //use sim::EnergyDeposits by default. This is heavy and may be dropped
+        if(sedlist.size()>0)
+        {
+          BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sedlist, fCaloPlane);
+        }
+        else //use sim::Channel -> IDE otherwise. This is usually kept but results in strange bugs.
+        {
+          BlipUtils::FillParticleInfo( *plist[i], pinfo[i], sIDElist, fCaloPlane);
+        }
         if( map_g4trkid_charge[pinfo[i].trackId] ) pinfo[i].numElectrons = (int)map_g4trkid_charge[pinfo[i].trackId];
+        else pinfo[i].numElectrons = pinfo[i].depElectrons; // JACOB ADDITION Jan22, 2026 for strange channel IDE results
         pinfo[i].index = i;
       }
       BlipUtils::MakeTrueBlips(pinfo, trueblips);
+      std::cout << "True blip size after make " << trueblips.size() << std::endl; 
       BlipUtils::MergeTrueBlips(trueblips, fTrueBlipMergeDist);
+      std::cout << "True blip size after merge " << trueblips.size() << std::endl; 
     }
 
-
     //=======================================
     // Map track IDs to the index in the vector
     //=======================================
@@ -647,7 +648,7 @@ namespace blip {
 
 
       // find associated track
-      if( fHitProducer == "gaushit" && fmtrk.isValid() ) {
+      if( fHitProducer == "specialblipgaushit" && fmtrk.isValid() ) {
         if(fmtrk.at(i).size()) hitinfo[i].trkid = fmtrk.at(i)[0]->ID();
 
       // if the hit collection didn't have associations made
@@ -664,8 +665,7 @@ namespace blip {
       if( hitinfo[i].trkid < 0 ) nhits_untracked++;
       //printf("  %lu   plane: %i,  wire: %i, time: %i\n",i,hitinfo[i].plane,hitinfo[i].wire,int(hitinfo[i].driftTime));
 
-    }//endloop over hits
-
+    }//endloop over hits    
     //for(auto& a : tpc_plane_hitsMap ) {
       //for(auto& b : a.second ) 
         //std::cout<<"TPC "<<a.first<<", plane "<<b.first<<": "<<b.second.size()<<" hits\n";
@@ -690,6 +690,7 @@ namespace blip {
 
 
     // Basic track inclusion cut: exclude hits that were tracked
+    int Tracked=0;
     for(size_t i=0; i<hitlist.size(); i++){
       if( hitinfo[i].trkid < 0 ) continue;
       auto it = map_trkid_index.find(hitinfo[i].trkid);
@@ -698,8 +699,10 @@ namespace blip {
       if( tracklist[trkindex]->Length() > fMaxHitTrkLength ) {
         hitIsTracked[i] = true;
         hitIsGood[i] = false;
+        Tracked++;
       }
     }
+    std::cout << Tracked << " hits were in tracks from pandora " << std::endl;
 
 
     // Filter based on hit properties. For hits that are a part of
@@ -731,7 +734,6 @@ namespace blip {
     // Hit clustering
     // ---------------------------------------------------
     std::map<int,std::map<int,std::vector<int>>> tpc_planeclustsMap;
-
     for(auto const& tpc_plane_hitsMap : cryo_tpc_plane_hitsMap ) {
 
     for(auto const& plane_hitsMap : tpc_plane_hitsMap.second ) {
@@ -892,7 +894,6 @@ namespace blip {
 
     float _matchQDiffLimit= (fMatchQDiffLimit <= 0 ) ? std::numeric_limits<float>::max() : fMatchQDiffLimit;
     float _matchMaxQRatio = (fMatchMaxQRatio  <= 0 ) ? std::numeric_limits<float>::max() : fMatchMaxQRatio;
-
     for(auto& tpcMap : tpc_planeclustsMap ) { // loop on TPCs
 
       //std::cout
@@ -1046,7 +1047,6 @@ namespace blip {
             for(auto& hc : hcGroup ) {
               hitclust[hc.ID].isMatched = true;
               for(auto hit : hitclust[hc.ID].HitIDs) hitinfo[hit].ismatch = true;
-
               // Diagnostic plots for successful 3-plane matches
               //if( picky && hc.Plane != fCaloPlane ) {
                 //float q1 = (float)newBlip.clusters[fCaloPlane].Charge;
@@ -1096,18 +1096,26 @@ namespace blip {
             // if we made it this far, the blip is good!
             // associate this blip with the hits and clusters within it
             newBlip.ID = blips.size();
-            blips.push_back(newBlip);
             for(auto& hc : hcGroup ) {
               hitclust[hc.ID].BlipID = newBlip.ID;
               for( auto& h : hc.HitIDs ) hitinfo[h].blipid = newBlip.ID;
             }
+            //BLIPS HAVE A COPY OF HITCLUSTERS NOT A POINTER
+            //UPDATE THE HITCLUSTER VARS THAT HAVE CHANGED SINCE CONSTRUCTION
+            for(int iclust=0; iclust<int(sizeof(newBlip.clusters)/sizeof(newBlip.clusters[0])); iclust++)
+            {
+              if(newBlip.clusters[iclust].ID>-1){
+                newBlip.clusters[iclust].BlipID = newBlip.ID;
+                newBlip.clusters[iclust].isMatched = true;
+              }
+            }
+            blips.push_back(newBlip);
 
 
           }//endif ncands > 0
         }//endloop over caloplane ("Plane A") clusters
       }//endif calo plane has clusters
     }//endloop over TPCs
-
     // Re-index the clusters after removing unmatched
     if( !keepAllClusts ) {
       std::vector<blip::HitClust> hitclust_filt;

sbndcode/BlipRecoSBND/BlipAna_module.cc

@@ -1582,6 +1582,9 @@ void BlipAna::PrintClusterInfo(const blip::HitClust& hc){
     hc.EdepID,
     hc.isMatched
   );
+  printf("G4 IDs contib \n");
+  for(int g4ID : hc.G4IDs) printf(" %i ", g4ID);
+  printf("\n");
 }
 
 void BlipAna::PrintBlipInfo(const blip::Blip& bl){

sbndcode/BlipRecoSBND/Utils/BlipUtils.cc

@@ -32,9 +32,8 @@ namespace BlipUtils {
   //===========================================================================
   // Provided a MCParticle, calculate everything we'll need for later calculations
   // and save into ParticleInfo object
-  void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){
-
-    // Get important info and do conversions
+void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo){
+  // Get important info and do conversions
     pinfo.particle    = part;
     pinfo.trackId     = part.TrackId();
     pinfo.isPrimary   = (int)(part.Process() == "primary");
@@ -50,17 +49,28 @@ namespace BlipUtils {
     pinfo.time        = /*ns ->mus*/1e-3 * part.T();
     pinfo.endtime     = /*ns ->mus*/1e-3 * part.EndT();
     pinfo.numTrajPts  = part.NumberTrajectoryPoints();
-
     // Pathlength (in AV) and start/end point
     pinfo.pathLength  = PathLength( part, pinfo.startPoint, pinfo.endPoint);
-
     // Central position of trajectory
     pinfo.position    = 0.5*(pinfo.startPoint+pinfo.endPoint);
-
     // Energy/charge deposited by this particle, found using SimEnergyDeposits 
     pinfo.depEnergy     = 0;
     pinfo.depElectrons  = 0;
+    return;
+}
+void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SEDVec_t& sedvec, int caloPlane){
+    FillParticleInfo( part, pinfo);
     for(auto& sed : sedvec ) {
+      if( -1*sed->TrackID() == part.TrackId() || sed->TrackID() == part.TrackId() ) {
+        pinfo.depEnergy     += sed->Energy();
+        pinfo.depElectrons  += sed->NumElectrons();
+      }
+    }
+    return;
+  }
+  void FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo, SIDEVec_t& sIDEvec, int caloPlane){
+    FillParticleInfo( part, pinfo);
+    for(auto& sed : sIDEvec ) {
       if( -1*sed.trackID == part.TrackId() || sed.trackID == part.TrackId() ) {
         pinfo.depEnergy     += sed.energy;
         pinfo.depElectrons  += sed.numElectrons;
@@ -103,8 +113,8 @@ namespace BlipUtils {
         for(size_t j=0; j<pinfo.size(); j++){
           simb::MCParticle& p = pinfo[j].particle;
           std::string pr = p.Process();
-          if( p.PdgCode() != 2112 && (pr == "eIoni" || pr == "muIoni" || pr == "hIoni") ){
-            if( IsAncestorOf(p.TrackId(),part.TrackId(),true) ) GrowTrueBlip(pinfo[j],tb);
+          if( p.PdgCode() != 2112 && p.PdgCode() != 22 && (pr == "eIoni" || pr == "muIoni" || pr == "hIoni") ){ //neutron and photons leave track
+            if( IsAncestorOf(p.TrackId(),part.TrackId(),true,true) ) GrowTrueBlip(pinfo[j],tb);
           }
         }
       }
@@ -151,7 +161,7 @@ namespace BlipUtils {
     // .. otherwise, check that the new particle
     // creation time is comparable to existing blip.
     // then calculate new energy-weighted position.
-    } else if ( fabs(tblip.Time-pinfo.time) < 3 ) {
+    } else if ( fabs(tblip.Time-pinfo.time) < 3) {
       float totE = tblip.Energy + pinfo.depEnergy;
       float w1 = tblip.Energy/totE;
       float w2 = pinfo.depEnergy/totE;
@@ -474,7 +484,7 @@ namespace BlipUtils {
   //====================================================================
   // Function to determine if a particle descended from another particle.
   // Allows option to break lineage at photons for contiguous parentage.
-  bool IsAncestorOf(int particleID, int ancestorID, bool breakAtPhots = false){
+  bool IsAncestorOf(int particleID, int ancestorID, bool breakAtPhots = false, bool breakAtNeutrons = false){
     art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
     const sim::ParticleList& plist = pi_serv->ParticleList();
     if( particleID == ancestorID  )       return true;
@@ -488,6 +498,7 @@ namespace BlipUtils {
       simb::MCParticle pM = pi_serv->TrackIdToParticle(p.Mother());
       if      ( pM.TrackId() == ancestorID )                      { return true;  }
       else if ( breakAtPhots == true && pM.PdgCode() == 22 )      { return false; }
+      else if ( breakAtNeutrons && pM.PdgCode() == 2112 )           { return false; }
       else if ( pM.Process() == "primary" || pM.TrackId() == 1 )  { return false; }
       else if ( pM.Mother() == 0 )                                { return false; }
       else    { particleID = pM.TrackId(); }              

sbndcode/BlipRecoSBND/Utils/BlipUtils.h

@@ -37,8 +37,8 @@
 #include "sbndcode/BlipRecoSBND/Utils/DataTypes.h"
 #include "TH1D.h"
 
-
-typedef std::vector<sim::IDE> SEDVec_t;
+typedef std::vector<art::Ptr<sim::SimEnergyDeposit>> SEDVec_t;
+typedef std::vector<sim::IDE> SIDEVec_t;
 
 geo::View_t kViews[3]={geo::kU, geo::kV, geo::kW};
 
@@ -48,7 +48,9 @@ namespace BlipUtils {
   // Functions related to blip reconstruction
   //###################################################
   //void      InitializeDetProps();
-  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane=2);
+  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SEDVec_t&, int plane);
+  void      FillParticleInfo(simb::MCParticle const&, blip::ParticleInfo&, SIDEVec_t&, int plane);
+  void      FillParticleInfo( const simb::MCParticle& part, blip::ParticleInfo& pinfo);
   //void      CalcPartDrift(blip::ParticleInfo&, int);
   //void      CalcTotalDep(float&,int&,float&, SEDVec_t&);
   void      MakeTrueBlips(std::vector<blip::ParticleInfo>&, std::vector<blip::TrueBlip>&); 
@@ -78,7 +80,7 @@ namespace BlipUtils {
   //bool    G4IdToMCTruth( int const, art::Ptr<simb::MCTruth>&);
   double  PathLength(const simb::MCParticle&, TVector3&, TVector3&);
   double  PathLength(const simb::MCParticle&);
-  bool    IsAncestorOf(int, int, bool);
+  bool    IsAncestorOf(int, int, bool, bool);
   double  DistToBoundary(const recob::Track::Point_t&);
   double  DistToLine(TVector3&, TVector3&, TVector3&);
   double  DistToLine2D(TVector2&, TVector2&, TVector2&);

sbndcode/BlipRecoSBND/blipreco_configs.fcl

@@ -8,7 +8,7 @@ sbnd_blipalg:
   HitProducer:          "specialblipgaushit"       #// input recob::Hits to use for blip reconstruction
   TrkProducer:          "blipPandoraTrackCopy"       #// input recob::Tracks to use for blip reconstruction
   GeantProducer:        "largeant"      #// input sim::MCParticles (getting true particle info)
-  SimEDepProducer:      "ionandscint"   #// input sim::SimEnergyDeposits (getting energy/electrons deposited)
+  SimEDepProducer:      "ionandscint:priorSCE:G4"   #// input sim::SimEnergyDeposits (getting energy/electrons deposited)
   SimChanProducer:      "simtpc2d:simpleSC:DetSim"      #// label for sim::SimChannels (getting drifted charge; optional)
   MaxHitTrkLength:      5.0             #// hits in trks > this length will be vetoed [cm]
   DoHitFiltering:       false           #// filter hits based on amp, width, RMS