Merge branch 'release/v10_12_02_01'

Author: nathanielerowe

CMakeLists.txt

@@ -15,7 +15,7 @@
 
 cmake_minimum_required(VERSION 3.20 FATAL_ERROR)
 
-set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.12.02)
+set(${PROJECT_NAME}_CMAKE_PROJECT_VERSION_STRING 10.12.02.01)
 find_package(cetmodules REQUIRED)
 project(sbndcode LANGUAGES CXX)
 

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:          "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)
-  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

sbndcode/CRT/CRTAna/CRTAnalysis_module.cc

@@ -63,6 +63,8 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
   // Required functions.
   void analyze(art::Event const& e) override;
 
+  void beginRun(art::Run const& r) override;
+
   void AnalysePTBs(std::vector<art::Ptr<raw::ptb::sbndptb>> &PTBVec);
 
   void AnalyseTDCs(std::vector<art::Ptr<sbnd::timing::DAQTimestamp>> &TDCVec);
@@ -95,8 +97,10 @@ class sbnd::crt::CRTAnalysis : public art::EDAnalyzer {
 
 private:
 
-  art::ServiceHandle<CRTGeoService> fCRTGeoService;
-  TPCGeoAlg fTPCGeoAlg;
+  art::ServiceHandle<CRTGeoService>              fCRTGeoService;
+  art::ServiceHandle<SBND::CRTChannelMapService> fCRTChannelMapService;
+
+  TPCGeoAlg         fTPCGeoAlg;
   CRTBackTrackerAlg fCRTBackTrackerAlg;
 
   std::string fMCParticleModuleLabel, fSimDepositModuleLabel, fFEBDataModuleLabel, fCRTStripHitModuleLabel,
@@ -640,37 +644,39 @@ sbnd::crt::CRTAnalysis::CRTAnalysis(fhicl::ParameterSet const& p)
       fTree->Branch("tdc_offset", "std::vector<uint64_t>", &_tdc_offset);
       fTree->Branch("tdc_name", "std::vector<std::string>", &_tdc_name);
     }
+}
 
+void sbnd::crt::CRTAnalysis::beginRun(art::Run const& r)
+{
   if(fDebug)
     {
-      for(auto const &[name, tagger] : fCRTGeoService->GetTaggers())
-        {
-          std::cout << "Tagger:  " << tagger.name << '\n'
-                    << "X - Min: " << tagger.minX << " Max: " << tagger.maxX << '\n'
-                    << "Y - Min: " << tagger.minY << " Max: " << tagger.maxY << '\n'
-                    << "Z - Min: " << tagger.minZ << " Max: " << tagger.maxZ << '\n' << std::endl;
-        }
+      std::cout << "\n==================================="
+                << "\nTaggers!"
+                << "\n===================================";
 
-      std::cout << std::endl;
+      for (auto const &[name, tagger] : fCRTGeoService->GetTaggers())
+        std::cout << std::endl << tagger;
 
-      for(auto const &[name, module] : fCRTGeoService->GetModules())
-        {
-          std::cout << "Module:  " << module.name << " (" << module.taggerName << ")" << '\n';
-          if(module.minos)
-            std::cout << "MINOS module" << std::endl;
-          std::cout << "X - Min: " << module.minX << " Max: " << module.maxX << " Diff: " << module.maxX - module.minX << '\n' 
-                    << "Y - Min: " << module.minY << " Max: " << module.maxY << " Diff: " << module.maxY - module.minY << '\n' 
-                    << "Z - Min: " << module.minZ << " Max: " << module.maxZ << " Diff: " << module.maxZ - module.minZ << '\n' 
-                    << "Orientation: " << module.orientation << '\n' << std::endl;
-        }
+      std::cout << "\n==================================="
+                << "\nModules!"
+                << "\n===================================";
 
-      std::cout << std::endl;
+      for (auto const &[name, module] : fCRTGeoService->GetModules())
+        std::cout << std::endl << module;
 
-      for(auto const &[name, sipm] : fCRTGeoService->GetSiPMs())
-        {
-          std::cout << "SiPM:  " << sipm.channel << " (" << sipm.channel/32 << " - " << sipm.channel%32 << ")" << '\n'
-                    << "x: " << sipm.x << " y: " << sipm.y << " z: " << sipm.z << std::endl;
-        }
+      std::cout << "\n==================================="
+                << "\nStrips!"
+                << "\n===================================";
+
+      for (auto const &[name, strip] : fCRTGeoService->GetStrips())
+        std::cout << std::endl << strip;
+
+      std::cout << "\n==================================="
+                << "\nSiPMs!"
+                << "\n===================================";
+
+      for (auto const &[name, sipm] : fCRTGeoService->GetSiPMs())
+        std::cout << std::endl << sipm;
     }
 }
 
@@ -1269,20 +1275,12 @@ void sbnd::crt::CRTAnalysis::AnalyseCRTClusters(const art::Event &e, const std::
       for(unsigned ii = 0; ii < _cl_nhits[i]; ++ii)
         {
           const auto striphit = striphits[ii];
-          _cl_channel_set[i][ii] = striphit->Channel();
-          _cl_adc_set[i][2*ii]   = striphit->ADC1();
-          _cl_adc_set[i][2*ii+1] = striphit->ADC2();
-          _cl_sh_ts0_set[i][ii]  = striphit->Ts0();
-          _cl_sh_ts1_set[i][ii]  = striphit->Ts1();
-
-          if(fDataMode)
-            {
-              art::ServiceHandle<SBND::CRTChannelMapService> ChannelMapService;
-              SBND::CRTChannelMapService::ModuleInfo_t module_info = ChannelMapService->GetModuleInfoFromOfflineID( striphit->Channel() / 32 );
-              _cl_sh_feb_mac5_set[i][ii] = ( module_info.valid ) ? module_info.feb_mac5 : 0;
-            }
-          else
-            _cl_sh_feb_mac5_set[i][ii] = striphit->Channel() / 32;
+          _cl_channel_set[i][ii]     = striphit->Channel();
+          _cl_adc_set[i][2*ii]       = striphit->ADC1();
+          _cl_adc_set[i][2*ii+1]     = striphit->ADC2();
+          _cl_sh_ts0_set[i][ii]      = striphit->Ts0();
+          _cl_sh_ts1_set[i][ii]      = striphit->Ts1();
+          _cl_sh_feb_mac5_set[i][ii] = fCRTChannelMapService->GetMAC5FromOfflineChannelID(striphit->Channel());
 
           /*
            * The below segment reimplements the CorrectTime() method

sbndcode/CRT/CRTAna/run_crtana.fcl

@@ -1,4 +1,4 @@
-#include "crtgeoservice_sbnd.fcl"
+#include "crt_services_sbnd.fcl"
 #include "crtana_sbnd.fcl"
 #include "simulationservices_sbnd.fcl"
 
@@ -8,12 +8,12 @@ services:
 {
   TFileService: { fileName: "crtana_sbnd.root" }
   @table::sbnd_basic_services
+  @table::crt_services_sbnd
   ParticleInventoryService:  @local::sbnd_particleinventoryservice
   BackTrackerService:        @local::sbnd_backtrackerservice
   DetectorClocksService:     @local::sbnd_detectorclocks
   LArPropertiesService:      @local::sbnd_properties
   DetectorPropertiesService: @local::sbnd_detproperties
-  CRTGeoService:             @local::crtgeoservice_sbnd
 }
 
 source: