Made associations with Space Points. Added cluster->Tagger() functionality. Added MC versus Data fcls

Author: aantonakis

sbndcode/CRT/CRTReco/CMakeLists.txt

@@ -34,5 +34,4 @@ simple_plugin(
               Eigen3::Eigen
 )
 
-
 install_fhicl()

sbndcode/CRT/CRTVeto/CMakeLists.txt

@@ -1,3 +1,5 @@
+
+
 art_make_library(
                  LIBRARIES
                         sbnobj::SBND_CRT

sbndcode/CRT/CRTVeto/CRTVetoProducer_module.cc

@@ -22,14 +22,15 @@
 
 #include "lardata/Utilities/AssociationUtil.h"
 
+#include "sbnobj/SBND/CRT/CRTEnums.hh"
 #include "sbnobj/SBND/CRT/CRTCluster.hh"
 #include "sbnobj/SBND/CRT/CRTSpacePoint.hh"
 #include "sbnobj/SBND/CRT/CRTTrack.hh"
 #include "sbnobj/SBND/CRT/CRTVeto.hh"
 
 #include "sbndcode/Geometry/GeometryWrappers/CRTGeoAlg.h"
 #include "sbndcode/CRT/CRTUtils/CRTCommonUtils.h"
-
+#include "lardata/Utilities/AssociationUtil.h"
 
 #include "TMath.h"
 
@@ -53,19 +54,24 @@ class sbnd::crt::CRTVetoProducer : public art::EDProducer {
   void produce(art::Event& e) override;
 
   // Functions to check if a point is in a specific Tagger
-  bool inSouth(const double &x, const double &y, const double &z);
-  bool inNorth(const double &x, const double &y, const double &z);
-  bool inEast(const double &x, const double &y, const double &z);
-  bool inWest(const double &x, const double &y, const double &z);
-  bool inTopLow(const double &x, const double &y, const double &z);
-  bool inTopHigh(const double &x, const double &y, const double &z);
+  bool inSouth(const CRTTagger t);
+  bool inNorth(const CRTTagger t);
+  bool inEast(const CRTTagger t);
+  bool inWest(const CRTTagger t);
+  bool inTopLow(const CRTTagger t);
+  bool inTopHigh(const CRTTagger t);
 
 private:
 
   CRTGeoAlg        fCRTGeoAlg;
   double           fWindowStart;
   double           fWindowEnd;
+  std::string      fCRTClusterModuleLabel;
   std::string      fCRTSpacePointModuleLabel;
+  bool		   fIsData;
+
+  std::vector<int> fVetoSpacePointIndices;
+  std::vector<art::Ptr<CRTSpacePoint>> fVetoSpacePoints;
 };
 
 
@@ -74,16 +80,35 @@ sbnd::crt::CRTVetoProducer::CRTVetoProducer(fhicl::ParameterSet const& p)
   , fCRTGeoAlg(p.get<fhicl::ParameterSet>("CRTGeoAlg"))
   , fWindowStart(p.get<double>("WindowStart"))
   , fWindowEnd(p.get<double>("WindowEnd"))
+  , fCRTClusterModuleLabel(p.get<std::string>("CRTClusterModuleLabel"))
   , fCRTSpacePointModuleLabel(p.get<std::string>("CRTSpacePointModuleLabel"))
+  , fIsData(p.get<bool>("IsData"))
 {
-  produces<CRTVeto>();
- // produces<art::Assns<CRTSpacePoint, CRTVeto>>();
+  produces<std::vector<CRTVeto>>();
+  produces<art::Assns<CRTVeto, CRTSpacePoint>>();
 }
 
 void sbnd::crt::CRTVetoProducer::produce(art::Event& e)
 {
-  //auto crtveto            = std::make_unique<CRTVeto>();
-  //auto vetoSpacePointAssn = std::make_unique<art::Assns<CRTSpacePoint, CRTVeto>>();
+
+  fVetoSpacePoints.clear();
+  fVetoSpacePointIndices.clear();
+ 
+  if (fIsData) {
+    std::cout << "Filling CRT Veto for Data" << std::endl;
+
+  }
+
+ 
+  // Get the Clusters for this event
+  art::Handle<std::vector<CRTCluster>> CRTClusterHandle;
+  e.getByLabel(fCRTClusterModuleLabel, CRTClusterHandle);
+  if(!CRTClusterHandle.isValid()){
+    std::cout << "CRTCluster product " << fCRTClusterModuleLabel << " not found..." << std::endl;
+    throw std::exception();
+  } 
+  std::vector<art::Ptr<CRTCluster>> CRTClusterVec;
+  art::fill_ptr_vector(CRTClusterVec, CRTClusterHandle);
 
   // Get the Space Points for this event
   art::Handle<std::vector<CRTSpacePoint>> CRTSpacePointHandle;
@@ -95,65 +120,77 @@ void sbnd::crt::CRTVetoProducer::produce(art::Event& e)
   std::vector<art::Ptr<CRTSpacePoint>> CRTSpacePointVec;
   art::fill_ptr_vector(CRTSpacePointVec, CRTSpacePointHandle);
 
+  art::FindManyP<CRTSpacePoint> cluster_sps(CRTClusterVec, e, fCRTSpacePointModuleLabel);
+  
+  // Get indices of Veto Space Points
+  /*
+  for (unsigned int i = 0; i < CRTSpacePointVec.size(); ++i) {
+
+    double t = CRTSpacePointVec.at(i)->Ts0()/1000; // Convert to us
+    if ((t > fWindowEnd) || (t < fWindowStart)) {
+      continue;
+    }
+    fVetoSpacePointIndices.push_back(i);
+
+  }
+  */
+
+  // initialize counters for veto hits in different taggers
   int S = 0;
   int N = 0;
   int E = 0;
   int W = 0;
   int TL = 0;
   int TH = 0;
 
-  // loop over the space points
-  for (unsigned int i = 0; i < CRTSpacePointVec.size(); i++){
-
-    double x = CRTSpacePointVec[i]->X();
-    double y = CRTSpacePointVec[i]->Y();
-    double z = CRTSpacePointVec[i]->Z();
-    double t = CRTSpacePointVec[i]->Ts0()/1000; // Convert to us
-    
-    // Debugging statement 
-    
-    /*
-    if ( (t <= 2) && (t >= 0.4) ) {
-      std::cout << std::endl;
-      std::cout << std::endl;
-      std::cout << "Found a CRT Hit within the South Wall Top Hat Window!" << std::endl;
+  // loop over clusters in this event
+  // Note: there are usually more clusters than space points
+  for(auto const &cluster : CRTClusterVec) {
+  
+    std::vector<art::Ptr<CRTSpacePoint>> clusterSpacePoints = cluster_sps.at(cluster.key());
+  
+    // loop over the space points in this cluster --> should be one per cluster
+    for (auto const &sp : clusterSpacePoints) {
+     
+      // Get the time and tagger for this space point
+      double t = sp->Ts0()/1000; // Convert to us
+      CRTTagger tag = cluster->Tagger();  
+
+      // check if in the beam window --> Make configureable
+      if ((t > fWindowEnd) || (t < fWindowStart)) {
+        continue;
+      }
       std::cout << std::endl;
+      std::cout << "Found Space Point in time with the beam !!!" << std::endl;
       std::cout << std::endl;
-
-    } 
-    */
-
-    // check if in the beam window --> Make configureable
-    if ((t > fWindowEnd) || (t < fWindowStart)) {
-      continue;
-    }
-    // We have a valid space point --> make assn
-    //util::CreateAssn(*this, e, *crtveto, CRTSpacePointVec[i], *vetoSpacePointAssn);
 
-    if (inSouth(x, y, z)) {
-      S += 1;
-    }
-    else if (inNorth(x, y, z)) {
-      N += 1;
-    }
-    else if (inEast(x, y, z)) {
-      E += 1;
-    }
-    else if (inWest(x, y, z)) {
-      W += 1;
-    }
-    else if (inTopLow(x, y, z)) {
-      TL += 1;
-    }
-    else if (inTopHigh(x, y, z)) {
-      TH += 1;
-    }
-    else {
-      std::cout << "Not a valid SpacePoint --> Outside Geometry !!!" << std::endl;
-    }
-
-  } // end loop over space points
-
+      fVetoSpacePoints.push_back(sp);
+      fVetoSpacePointIndices.push_back(sp.key());
+
+      if (inSouth(tag)) {
+        S += 1;
+      }
+      else if (inNorth(tag)) {
+        N += 1;
+      }
+      else if (inEast(tag)) {
+        E += 1;
+      }
+      else if (inWest(tag)) {
+        W += 1;
+      }
+      else if (inTopLow(tag)) {
+        TL += 1;
+      }
+      else if (inTopHigh(tag)) {
+        TH += 1;
+      }
+      else {
+        std::cout << "Not a valid SpacePoint --> Outside Geometry !!!" << std::endl;
+      }
+    } // end loop over space points
+  } // end of loop over clusters
+ 
   int Nwalls_all = S + N + E + W;
   int Nwalls_noNorth = S + E + W;
   int Ntop = TL + TH;
@@ -185,58 +222,87 @@ void sbnd::crt::CRTVetoProducer::produce(art::Event& e)
   }
 
   //crtveto =  CRTVeto(v0, v1, v2, v3); 
+  //
+  //
+  // Make Products for the Event
   auto crtveto            = std::make_unique<CRTVeto>(v0, v1, v2, v3, v4);
- 
-  e.put(std::move(crtveto));
-  //e.put(std::move(vetoSpacePointAssn));
+  auto vetoVec            = std::make_unique<std::vector<CRTVeto>>();
+  vetoVec->push_back(*crtveto);
+
+  auto vetoSpacePointAssn = std::make_unique<art::Assns<CRTVeto, CRTSpacePoint>>();
+
+  //auto const vetoPtrMaker = art::PtrMaker<CRTVeto>(e);
+  //auto const vetoPtr = vetoPtrMaker(0); 
+  //auto const CRTSpacePointPtrMaker = art::PtrMaker<CRTSpacePoint>(e, CRTSpacePointHandle.id());
+
+  if (fVetoSpacePoints.size() != fVetoSpacePointIndices.size()) {
+    std::cout << std::endl;
+    std::cout << "Problem: Number of Space Points does not match!" << std::endl;
+    std::cout << std::endl;
+  }
+
+  util::CreateAssn(*this, e, *vetoVec, fVetoSpacePoints, *vetoSpacePointAssn);
+  
+  /*
+  for (unsigned int i = 0; i < fVetoSpacePoints.size(); ++i) {
+    std::cout << "Veto SP index " << fVetoSpacePointIndices.at(i) << std::endl;
+    //auto const CRTSpacePointPtr = CRTSpacePointPtrMaker(fVetoSpacePointIndices.at(i));
+    //vetoSpacePointAssn->addSingle(vetoPtr, CRTSpacePointPtr);
+    //vetoSpacePointAssn->addSingle(crtveto, CRTSpacePointPtr);
+    util::CreateAssn(*this, e, *vetoVec, fVetoSpacePoints.at(i), *vetoSpacePointAssn);
+  }
+  */
+
+  //e.put(std::move(crtveto));
+  e.put(std::move(vetoVec));
+  e.put(std::move(vetoSpacePointAssn));
 
 } // end of produce
 
-bool sbnd::crt::CRTVetoProducer::inSouth(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inSouth(const CRTTagger t)
 {
-  if ((z > -1000) && (z < -179) && (y > -370) && (y < 400)) {
+  if (t == kSouthTagger) {
     return true;
   }
   return false;
 }
-bool sbnd::crt::CRTVetoProducer::inNorth(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inNorth(const CRTTagger t)
 {
-  // check that this isn't affected by x from east and west --> TODO
-  if ((z > 746) && (y > -370) && (y < 400)) {
+  if (t == kNorthTagger) {
     return true;
   }
   return false;
 }
 
-bool sbnd::crt::CRTVetoProducer::inEast(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inEast(const CRTTagger t)
 {
   //east_sel = "x < -380 and z < 770 and  -370 < y < 400"
-  if ((x < -380) && (z < 770) && (y > -370) && (y < 400)) {
+  if (t == kEastTagger) {
     return true;
   }
   return false;
 }
 
-bool sbnd::crt::CRTVetoProducer::inWest(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inWest(const CRTTagger t)
 {
   //west_sel = "x < -380 and z < 770 and  -370 < y < 400"
-  if ((x > 380) && (z < 770) && (y > -370) && (y < 400)) {
+  if (t == kWestTagger) {
     return true;
   }
   return false;
 }
 
-bool sbnd::crt::CRTVetoProducer::inTopLow(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inTopLow(const CRTTagger t)
 {
-  if ((y > 400) && (y < 700)) {
+  if (t == kTopLowTagger) {
     return true;
   }
   return false;
 }
 
-bool sbnd::crt::CRTVetoProducer::inTopHigh(const double &x, const double &y, const double &z)
+bool sbnd::crt::CRTVetoProducer::inTopHigh(const CRTTagger t)
 {
-  if (y > 700) {
+  if (t == kTopHighTagger) {
     return true;
   }
   return false;

sbndcode/CRT/CRTVeto/crtvetoproducer_sbnd.fcl

@@ -7,14 +7,19 @@ BEGIN_PROLOG
 crtvetoproducer_sbnd:
 {
    CRTGeoAlg:                  @local::crtgeoalg_sbnd
+   CRTClusterModuleLabel:      "crtclustering"
    CRTSpacePointModuleLabel:   "crtspacepoints"
-   WindowStart:                0.4
-   WindowEnd:                  2.0
+   WindowStart:                0.265
+   WindowEnd:                  1.965
+   IsData:                     false
    module_type:                "CRTVetoProducer"
 }
 
 crtvetoproducer_data_sbnd: @local::crtvetoproducer_sbnd
 crtvetoproducer_data_sbnd.CRTGeoAlg:     @local::crtgeoalg_data_sbnd
+crtvetoproducer_data_sbnd.IsData:        true
+crtvetoproducer_data_sbnd.WindowStart:   0.365
+crtvetoproducer_data_sbnd.WindowEnd:     2.065
 
 
 END_PROLOG

sbndcode/CRT/CRTVeto/run_crtveto_data.fcl

@@ -0,0 +1,10 @@
+#include "crt_channel_map_service.fcl"
+#include "crt_calib_service.fcl"
+#include "run_crtveto.fcl"
+
+services.message:              @local::sbnd_message_services
+services.CRTChannelMapService: @local::crt_channel_map_standard
+services.CRTCalibService:      @local::crt_calib_service
+
+physics.producers.crtveto:     @local::crtvetoproducer_data_sbnd
+