Blip True EDep from sim:IDE Object

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc

@@ -357,18 +357,32 @@ namespace blip {
     if (evt.getByLabel(fGeantProducer,pHandle))
       art::fill_ptr_vector(plist, pHandle);
 
-    // -- 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);
-    }
-
+    // -- 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)
     art::Handle<std::vector<sim::SimChannel> > simchanHandle;
     std::vector<art::Ptr<sim::SimChannel> > simchanlist;
-    if (evt.getByLabel(fSimChanProducer,simchanHandle)) 
+    if (evt.getByLabel(fSimChanProducer,simchanHandle))
+      { 
       art::fill_ptr_vector(simchanlist, simchanHandle);
+      //Loop over channels to get full sedlist
+      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 &
+	    }
+	}
+      }
 
     // -- hits (from input module, usually track-masked subset of gaushit)
     art::Handle< std::vector<recob::Hit> > hitHandle;

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h

@@ -36,6 +36,7 @@
 #include "lardataobj/AnalysisBase/BackTrackerMatchingData.h"
 #include "larsim/MCCheater/ParticleInventoryService.h"
 #include "larcore/Geometry/Geometry.h"
+#include "larcorealg/Geometry/WireReadoutGeom.h"
 #include "larcorealg/Geometry/GeometryCore.h"
 #include "larreco/Calorimetry/CalorimetryAlg.h"
 #include "art/Framework/Principal/Event.h"
@@ -117,6 +118,7 @@ namespace blip {
     float   kDriftVelocity;
     float   kTickPeriod;
     int     kNumChannels;
+    int     fCaloPlane;
 
    private:
 
@@ -172,7 +174,7 @@ namespace blip {
     bool                fKeepAllClusts[kNplanes];
 
     // --- Calorimetry configs ---
-    int                 fCaloPlane;
+    //int                 fCaloPlane;
     float               fCalodEdx;
     float               fESTAR_p0;
     float     	      	fESTAR_p1;

sbndcode/BlipRecoSBND/BlipRecoProducer_module.cc

@@ -88,6 +88,7 @@ BlipRecoProducer::BlipRecoProducer(fhicl::ParameterSet const & pset)
   produces< std::vector<  recob::SpacePoint > >();
   produces< art::Assns <  recob::Hit, recob::SpacePoint> >();
   produces< std::vector< blip::Blip > >(); 
+  produces< std::vector<blip::HitClust> >();
 
   //produces< art::Assns <  blip::Blip, recob::SpacePoint > >();
   produces< art::Assns <  blip::Blip,         recob::Hit> >();
@@ -124,6 +125,7 @@ void BlipRecoProducer::produce(art::Event & evt)
   std::unique_ptr< art::Assns <blip::Blip, recob::Hit> >  assn_blip_hit_v(std::make_unique<art::Assns<blip::Blip, recob::Hit> >() );
   std::unique_ptr< std::vector< recob::SpacePoint> > SpacePoint_v(std::make_unique<std::vector<recob::SpacePoint>>());
   std::unique_ptr< art::Assns <recob::Hit, recob::SpacePoint> >  assn_hit_sps_v(std::make_unique<art::Assns<recob::Hit,recob::SpacePoint>>() );
+  std::unique_ptr< std::vector< blip::HitClust> > collection_hitclust(std::make_unique<std::vector<blip::HitClust>>());
 
 
   art::PtrMaker<blip::Blip> makeBlipPtr(evt);
@@ -180,6 +182,12 @@ void BlipRecoProducer::produce(art::Event & evt)
     }
 
   }
+  //adding all the collection hit clusters
+  for(int iclust=0; iclust<int(fBlipAlg->hitclust.size()); iclust++)
+  {
+    if( (fBlipAlg->hitclust)[iclust].Plane != fBlipAlg->fCaloPlane) continue;
+    collection_hitclust->push_back((fBlipAlg->hitclust)[iclust]);
+  }
 
   //===========================================
   // Put them on the event
@@ -189,6 +197,7 @@ void BlipRecoProducer::produce(art::Event & evt)
   evt.put(std::move(blip_v));
   //evt.put(std::move(assn_blip_sps_v));
   evt.put(std::move(assn_blip_hit_v));
+  evt.put(std::move(collection_hitclust));
 }//END EVENT LOOP
 
 DEFINE_ART_MODULE(BlipRecoProducer)

sbndcode/BlipRecoSBND/Utils/BlipUtils.cc

@@ -61,12 +61,11 @@ namespace BlipUtils {
     pinfo.depEnergy     = 0;
     pinfo.depElectrons  = 0;
     for(auto& sed : sedvec ) {
-      if( sed->TrackID() == part.TrackId() ) {
-        pinfo.depEnergy     += sed->Energy();
-        pinfo.depElectrons  += sed->NumElectrons();
+      if( -1*sed.trackID == part.TrackId() || sed.trackID == part.TrackId() ) {
+        pinfo.depEnergy     += sed.energy;
+        pinfo.depElectrons  += sed.numElectrons;
       }
     }
-
     return;
 
   }
@@ -86,7 +85,7 @@ namespace BlipUtils {
       //std::cout<<"Making true blip for "<<part.TrackId()<<" (PDG "<<part.PdgCode()<<", which deposited "<<pinfo[i].depEnergy<<"\n";
 
       // If this is a photon or neutron, don't even bother!
-      if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) continue;
+      //if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) continue;
 
       // If this is an electron that came from another electron, it would 
       // have already been grouped as part of the contiguous "blip" previously.
@@ -139,7 +138,7 @@ namespace BlipUtils {
     simb::MCParticle& part = pinfo.particle;
 
     // Skip neutrons, photons
-    if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) return;
+    // if( part.PdgCode() == 2112 || part.PdgCode() == 22 ) return;
 
     // Check that path length isn't zero
     if( !pinfo.pathLength ) return;

sbndcode/BlipRecoSBND/Utils/BlipUtils.h

@@ -38,7 +38,7 @@
 #include "TH1D.h"
 
 
-typedef std::vector<art::Ptr<sim::SimEnergyDeposit>> SEDVec_t;
+typedef std::vector<sim::IDE> SEDVec_t;
 
 geo::View_t kViews[3]={geo::kU, geo::kV, geo::kW};
 

sbndcode/BlipRecoSBND/Utils/classes_def.xml

@@ -3,6 +3,8 @@
 <class name="blip::Blip"/>
 <class name="std::vector<blip::Blip>"/>
 <class name="blip::HitClust"/>
+<class name="std::vector<blip::HitClust>"/>
+<class name="art::Wrapper<std::vector<blip::HitClust> >"/>
 <class name="blip::TrueBlip"/>
 <class name="std::map<int,TVector3>"/>
 <class name="art::Assns<blip::Blip,recob::Hit,void>"/>

sbndcode/BlipRecoSBND/blipreco_configs.fcl

@@ -9,7 +9,7 @@ sbnd_blipalg:
   TrkProducer:          "pandoraTrack"       #// 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)
-  SimChanProducer:      "simdrift"      #// label for sim::SimChannels (getting drifted charge; optional)
+  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
   HitClustWireRange:    1               #// clustering wire width