Feature/adding blip to caf

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.cc

@@ -1072,7 +1072,7 @@ namespace blip {
               TVector3 p2(b.X(), b.Y(), b.Z() );
               // TO-DO: if this track starts or ends at a TPC boundary, 
               // we should extend p1 or p2 to outside the AV to avoid blind spots
-              TVector3 bp = newBlip.Position;
+              TVector3 bp(newBlip.Position.X(), newBlip.Position.Y(), newBlip.Position.Z());
               float d = BlipUtils::DistToLine(p1,p2,bp);
               if( d > 0 ) {
                 // update closest trkdist

sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h

@@ -40,7 +40,8 @@
 #include "larcorealg/Geometry/GeometryCore.h"
 #include "larreco/Calorimetry/CalorimetryAlg.h"
 #include "art/Framework/Principal/Event.h"
-
+#include "larcore/Geometry/WireReadout.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
 
 // Microboone includes
 //#include "ubevt/Database/TPCEnergyCalib/TPCEnergyCalibService.h"
@@ -54,6 +55,7 @@
 
 // Blip-specific utils
 #include "sbndcode/BlipRecoSBND/Utils/BlipUtils.h"
+#include "sbnobj/SBND/Blip/BlipDataTypes.h"
 
 // ROOT stuff
 #include "TH1D.h"

sbndcode/BlipRecoSBND/Alg/CMakeLists.txt

@@ -11,6 +11,7 @@ art_make_library(LIBRARY_NAME  sbndcode_BlipRecoAlg
                            lardataobj::RawData
                            lardataobj::RecoBase
                            lardata::RecoObjects
+			               larreco::Calorimetry
                            larpandora::LArPandoraInterface
                            nusimdata::SimulationBase
                            cetlib::cetlib
@@ -33,8 +34,8 @@ art_make_library(LIBRARY_NAME  sbndcode_BlipRecoAlg
                            ROOT::Gdml
                            sbndcode_CRTUtils
                            sbnobj::Common_CRT
-                           #sbndcode_CosmicIdUtils
                            sbndcode_BlipUtils
+			               sbnobj::SBND_Blip
 )
 
 install_headers()

sbndcode/BlipRecoSBND/BlipAna_module.cc

@@ -32,6 +32,7 @@
 #include "lardataobj/AnalysisBase/Calorimetry.h"
 #include "larevt/SpaceChargeServices/SpaceChargeService.h"
 #include "cetlib/search_path.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
 
 // SBND-specific includes
 #include "sbndcode/BlipRecoSBND/Alg/BlipRecoAlg.h"
@@ -75,7 +76,8 @@ const int kMaxTrks    =   1000;
 const int kMaxBlips   =   5000;
 const int kMaxG4      =   30000;
 const int kMaxEDeps   =   10000;
-const int kMaxTrkPts  =   2000;  
+const int kMaxTrkPts  =   2000; 
+const int kNplanes = blip::kNplanes;
 
 class BlipAna;
 

sbndcode/BlipRecoSBND/BlipRecoProducer_module.cc

@@ -33,6 +33,7 @@
 #include "lardataobj/RawData/RawDigit.h"
 #include "lardataobj/RawData/raw.h"
 #include "lardata/Utilities/AssociationUtil.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
 
 // C++ includes
 #include <vector>

sbndcode/BlipRecoSBND/CMakeLists.txt

@@ -8,7 +8,7 @@ set(      MODULE_LIBRARIES larcorealg::Geometry
                            larevt::Filters
                            lardataobj::RawData
                            lardataobj::RecoBase
-			                     larreco::Calorimetry
+			               larreco::Calorimetry
                            larreco::RecoAlg
                            lardata::RecoObjects
                            larpandora::LArPandoraInterface
@@ -34,12 +34,12 @@ set(      MODULE_LIBRARIES larcorealg::Geometry
                            sbndcode_RecoUtils
                            sbndcode_OpDetSim
                            sbndcode_BlipUtils
+                           sbnobj::SBND_Blip
                            sbndcode_BlipRecoAlg
 )
 
-add_subdirectory(Utils)
 add_subdirectory(Alg)
-
+add_subdirectory(Utils)
 cet_build_plugin(BlipAna art::Module SOURCE BlipAna_module.cc LIBRARIES ${MODULE_LIBRARIES} )
 cet_build_plugin(BlipRecoProducer art::Module SOURCE BlipRecoProducer_module.cc LIBRARIES ${MODULE_LIBRARIES} )
 

sbndcode/BlipRecoSBND/Utils/BlipUtils.cc

@@ -55,7 +55,9 @@ namespace BlipUtils {
     pinfo.pathLength  = PathLength( part, pinfo.startPoint, pinfo.endPoint);
 
     // Central position of trajectory
-    pinfo.position    = 0.5*(pinfo.startPoint+pinfo.endPoint);
+    pinfo.position.SetXYZ(0.5*(pinfo.startPoint.X()+pinfo.endPoint.X()),
+                         0.5*(pinfo.startPoint.Y()+pinfo.endPoint.Y()),
+                         0.5*(pinfo.startPoint.Z()+pinfo.endPoint.Z()) );
 
     // Energy/charge deposited by this particle, found using SimEnergyDeposits 
     pinfo.depEnergy     = 0;
@@ -155,7 +157,9 @@ namespace BlipUtils {
       float totE = tblip.Energy + pinfo.depEnergy;
       float w1 = tblip.Energy/totE;
       float w2 = pinfo.depEnergy/totE;
-      tblip.Position    = w1*tblip.Position + w2*pinfo.position;
+      tblip.Position.SetXYZ( w1*tblip.Position.X() + w2*pinfo.position.X(), 
+                            w1*tblip.Position.Y() + w2*pinfo.position.Y(),
+                            w1*tblip.Position.Z() + w2*pinfo.position.Z());
       tblip.Time        = w1*tblip.Time     + w2*pinfo.time;
       tblip.LeadCharge  = pinfo.depElectrons;
     // ... if the particle isn't a match, show's over
@@ -196,15 +200,17 @@ namespace BlipUtils {
         // check that the times are similar (we don't want to merge
         // together a blip that happened much later but in the same spot)
         if( fabs(blip_i.Time - blip_j.Time) > 5 ) continue;
-        float d = (blip_i.Position-blip_j.Position).Mag();
+        float d = TMath::Sqrt((blip_i.Position-blip_j.Position).Mag2());
         if( d < dmin ) {
           isGrouped.at(j) = true;
           //float totE = blip_i.Energy + blip_j.Energy;
           float totQ = blip_i.DepElectrons + blip_j.DepElectrons;
           float w1 = blip_i.DepElectrons/totQ;
           float w2 = blip_j.DepElectrons/totQ;
           blip_i.Energy       += blip_j.Energy;
-          blip_i.Position     = w1*blip_i.Position + w2*blip_j.Position;
+          blip_i.Position.SetXYZ( w1*blip_i.Position.X() + w2*blip_j.Position.X(), 
+                                  w1*blip_i.Position.Y() + w2*blip_j.Position.Y(),
+                                  w1*blip_i.Position.Z() + w2*blip_j.Position.Z());
           blip_i.DriftTime    = w1*blip_i.DriftTime+ w2*blip_j.DriftTime; 
           blip_i.Time         = w1*blip_i.Time + w2*blip_j.Time; 
           blip_i.DepElectrons += blip_j.DepElectrons;
@@ -419,12 +425,16 @@ namespace BlipUtils {
     // YZ-plane, as well as the mean difference between intersection points.
     newblip.Position.SetXYZ(0,0,0);
     if( wirex.size() == 1 ) {
-      newblip.Position= wirex[0];
+      newblip.Position = wirex[0];
     } else {
       newblip.SigmaYZ = 0;
       double fact = 1./wirex.size();
-      for(auto& v : wirex ) newblip.Position  += v * fact;
-      for(auto& v : wirex ) newblip.SigmaYZ   += (v-newblip.Position).Mag() * fact;
+      for(auto& v : wirex ) newblip.Position.SetXYZ( newblip.Position.X()  + v.X() * fact, 
+                                                     newblip.Position.Y()  + v.Y() * fact, 
+                                                     newblip.Position.Z()  + v.Z() * fact);
+      for(auto& v : wirex ) newblip.SigmaYZ   += TMath::Sqrt( pow(v.X()-newblip.Position.X(), 2) + 
+                                                              pow(v.Y()-newblip.Position.Y(), 2) + 
+                                                              pow(v.Z()-newblip.Position.Z(), 2)) * fact;
       // Ensure that difference between intersection points is
       // consistent with the maximal wire extent
       if( newblip.SigmaYZ > std::max(1.,0.5*newblip.dYZ) ) return newblip;
@@ -613,7 +623,7 @@ namespace BlipUtils {
 
   //=============================================================================
   // Length of particle trajectory
-  double PathLength(const simb::MCParticle& part, TVector3& start, TVector3& end)
+  double PathLength(const simb::MCParticle& part, geo::Point_t& start, geo::Point_t& end)
   {
     int n = part.NumberTrajectoryPoints();
     if( n <= 1 ) return 0.;
@@ -633,7 +643,7 @@ namespace BlipUtils {
     return L;
   }
   double PathLength(const simb::MCParticle& part){
-    TVector3 a,b;
+    geo::Point_t a,b;
     return PathLength(part,a,b);
   }
 

sbndcode/BlipRecoSBND/Utils/BlipUtils.h

@@ -27,14 +27,15 @@
 #include "larcore/CoreUtils/ServiceUtil.h"
 #include "larcore/Geometry/Geometry.h"
 #include "larcore/Geometry/WireReadout.h"
+#include "larcoreobj/SimpleTypesAndConstants/geo_vectors.h"
 //#include "larcorealg/Geometry/GeometryCore.h"
 //#include "larcorealg/Geometry/WireReadoutGeom.h"
 
 // c++
 #include <vector>
 #include <map>
 
-#include "sbndcode/BlipRecoSBND/Utils/DataTypes.h"
+#include "sbnobj/SBND/Blip/BlipDataTypes.h"
 #include "TH1D.h"
 
 
@@ -76,7 +77,7 @@ namespace BlipUtils {
   //void    HitTruth(art::Ptr<recob::Hit> const&, int&, float&, float&, float&);
   //si_t    HitTruthIds( art::Ptr<recob::Hit> const&);
   //bool    G4IdToMCTruth( int const, art::Ptr<simb::MCTruth>&);
-  double  PathLength(const simb::MCParticle&, TVector3&, TVector3&);
+  double  PathLength(const simb::MCParticle&, geo::Point_t&, geo::Point_t&);
   double  PathLength(const simb::MCParticle&);
   bool    IsAncestorOf(int, int, bool);
   double  DistToBoundary(const recob::Track::Point_t&);

sbndcode/BlipRecoSBND/Utils/CMakeLists.txt

@@ -4,7 +4,6 @@ cet_make_library(
                 LIBRARY_NAME  
                 sbndcode_BlipUtils
                 LIBRARIES 
-                PUBLIC
                 larcorealg::Geometry
                 larcore::Geometry_Geometry_service
                 larsim::Simulation
@@ -38,7 +37,7 @@ cet_make_library(
                 ROOT::Gdml
                 sbndcode_CRTUtils
                 sbnobj::Common_CRT
-                #sbndcode_CosmicIdUtils
+		sbnobj::SBND_Blip
 )
 
 art_dictionary(DICTIONARY_LIBRARIES sbndcode_BlipUtils)