yz norm feature

sbndcode/Calibration/CMakeLists.txt

@@ -2,3 +2,4 @@ add_subdirectory(OpReco)
 add_subdirectory(CRT)
 add_subdirectory(DQM)
 add_subdirectory(PDSDatabaseInterface)
+add_subdirectory(TPCCalorimetry)

sbndcode/Calibration/TPCCalorimetry/CMakeLists.txt

@@ -0,0 +1,22 @@
+
+set(TOOL_LIBRARIES
+  lardataobj::RecoBase
+  larreco::Calorimetry
+  larcorealg::Geometry
+  lardataalg::DetectorInfo
+  lardata::ArtDataHelper
+  canvas::canvas
+  art::Framework_Services_Registry
+  art::Utilities
+  fhiclcpp::fhiclcpp
+  cetlib::cetlib
+  cetlib_except::cetlib_except
+  ROOT::Core
+  ROOT::Hist
+)
+
+cet_build_plugin(NormalizeYZ art::tool LIBRARIES ${TOOL_LIBRARIES})
+
+install_headers()
+install_fhicl()
+install_source()

sbndcode/Calibration/TPCCalorimetry/NormalizeYZ_tool.cc

@@ -0,0 +1,143 @@
+
+// Framework Includes
+#include "art/Framework/Core/EDProducer.h"
+#include "art/Framework/Principal/Event.h"
+#include "art/Framework/Principal/Handle.h"
+#include "art/Framework/Services/Registry/ServiceHandle.h"
+#include "art/Utilities/ToolMacros.h"
+#include "cetlib/cpu_timer.h"
+#include "fhiclcpp/ParameterSet.h"
+#include "messagefacility/MessageLogger/MessageLogger.h"
+
+#include "canvas/Persistency/Common/Ptr.h"
+
+// Tool include
+#include "larreco/Calorimetry/INormalizeCharge.h"
+
+// Services
+#include "lardata/DetectorInfoServices/DetectorClocksService.h"
+
+// ROOT includes
+#include "TFile.h"
+#include "TH2F.h"
+
+// C++ includes
+#include <map>
+#include <string>
+#include <stdexcept>
+#include <iostream>
+
+
+
+namespace sbnd {
+  namespace calo {
+
+    class NormalizeYZ : public INormalizeCharge{
+    public:
+      explicit NormalizeYZ(fhicl::ParameterSet const& pset);
+
+      void configure(const fhicl::ParameterSet& pset) override;
+
+      double Normalize(double dQdx,
+                       const art::Event& evt,
+                       const recob::Hit& hit,
+                       const geo::Point_t& location,
+                       const geo::Vector_t& direction,
+                       double t0) override;
+
+    private:
+      void reconfigure(const fhicl::ParameterSet& pset);
+
+      std::map<std::string, TH2F*> fCorrHists;
+      std::string fFileName;
+      bool fVerbose;
+    };
+
+    NormalizeYZ::NormalizeYZ(fhicl::ParameterSet const& pset){
+      reconfigure(pset);  // delegate config to reusable function
+    }
+
+    void NormalizeYZ::configure(const fhicl::ParameterSet& pset){
+      reconfigure(pset); 
+    }
+
+    void NormalizeYZ::reconfigure(const fhicl::ParameterSet& pset){
+      fFileName = pset.get<std::string>("FileName");
+      fVerbose = pset.get<bool>("Verbose", false);
+
+      std::string fname;
+      cet::search_path sp("FW_SEARCH_PATH");
+      sp.find_file(fFileName, fname);  
+
+      TFile* f = TFile::Open(fname.c_str(), "READ");
+      if (!f || f->IsZombie()) {
+        throw cet::exception("NormalizeYZ") << "Failed to open correction map file: " << fFileName;
+      }
+
+      for (int plane = 0; plane < 3; plane++){   // planes : 2 inductions (idx : 0, 1) and 1 collection (idx : 2)
+        for (int tpc = 0; tpc < 2; tpc++){     // tpc : east (idx : 0) and west (idx : 1) TPCs
+          std::string histname = Form("CzyHist_%d_%d", plane, tpc);
+          TH2F* h = (TH2F*)f->Get(histname.c_str());
+          if (!h) {
+            throw cet::exception("NormalizeYZ") << "Missing histogram in file: " << histname;
+          }
+
+          fCorrHists[Form("plane%d_%d", plane, tpc)] = (TH2F*)h->Clone();
+          fCorrHists[Form("plane%d_%d", plane, tpc)]->SetDirectory(nullptr);
+        }
+      }
+      f->Close();
+    }
+
+    double NormalizeYZ::Normalize(double dQdx,
+				  const art::Event& evt,
+				  const recob::Hit& hit,
+				  const geo::Point_t& location,
+				  const geo::Vector_t&,
+				  double){
+
+      int plane = hit.WireID().Plane;
+      //int tpc = hit.WireID().TPC;
+
+      // just to be sure and consistent with the logic of input YZ map
+      // seems like current setup of directly calling hit.WireID().TPC has a tolerance of 30cm
+      // meaning - an approx region of -30<x<30 lies in both the TPCs    
+      int tpc;
+      if(location.X()<0){
+	tpc = 0;
+      } else tpc = 1;
+
+
+      std::string key = Form("plane%d_%d", plane, tpc);
+
+      auto it = fCorrHists.find(key);
+      if (it == fCorrHists.end()) {
+        mf::LogWarning("NormalizeYZ") << "No correction histogram for " << key << ". Returning uncorrected dQdx";
+        return dQdx;
+      }
+
+      TH2F* hCorr = it->second;
+      int binX = hCorr->GetXaxis()->FindBin(location.Z());
+      int binY = hCorr->GetYaxis()->FindBin(location.Y());
+      double scale = hCorr->GetBinContent(binX, binY);
+
+      if (fVerbose){
+        std::cout << "[NormalizeYZ] Plane: " << plane << ", TPC: " << tpc
+                  << ", Y: " << location.Y() << ", Z: " << location.Z()
+                  << ", Scale: " << scale << ", dQdx (raw): " << dQdx
+                  << ", dQdx (corrected): " << dQdx / scale << std::endl;
+      }
+
+      if (scale < 1e-3){
+        mf::LogWarning("NormalizeYZ") << "Invalid scale at (Y,Z)=(" << location.Y() << "," << location.Z() << "). Returning uncorrected dQdx";
+        return dQdx;
+      }
+
+      return dQdx / scale;
+    }
+
+    DEFINE_ART_CLASS_TOOL(NormalizeYZ)
+
+  } // namespace calo
+} // namespace sbnd
+

sbndcode/Calibration/TPCCalorimetry/normtools_sbnd.fcl

@@ -0,0 +1,19 @@
+BEGIN_PROLOG
+
+yznorm_hist_data: {
+  tool_type: NormalizeYZ
+  FileName: "YZmaps/yz_correction_map_data1e20.root"
+  Verbose: false
+}
+
+yznorm_hist_mc: {
+  tool_type: NormalizeYZ
+  FileName: "YZmaps/yz_correction_map_mcp2025b5e18.root"
+  Verbose: false
+}
+
+# list of normalization tools - by far only YZ correction is implemented
+sbnd_calonormtoolsdata: [@local::yznorm_hist_data]
+sbnd_calonormtoolsmc: [@local::yznorm_hist_mc]
+
+END_PROLOG

sbndcode/LArSoftConfigurations/calorimetry_sbnd.fcl

@@ -1,4 +1,5 @@
 #include "calorimetry.fcl"
+#include "normtools_sbnd.fcl"
 
 BEGIN_PROLOG
 
@@ -21,11 +22,13 @@ sbnd_gnewcalomc:              @local::standard_gnocchicalo
 sbnd_gnewcalomc.CaloAlg:      @local::sbnd_calorimetryalgmc
 sbnd_gnewcalomc.ChargeMethod: 3
 sbnd_gnewcalomc.SpacePointModuleLabel: @erase
+sbnd_gnewcalomc.NormTools: @local::sbnd_calonormtoolsmc   # YZ normalization for MC
 
 sbnd_gnewcalodata:              @local::standard_gnocchicalo
 sbnd_gnewcalodata.CaloAlg:      @local::sbnd_calorimetryalgdata
 sbnd_gnewcalodata.ChargeMethod: 3
 sbnd_gnewcalodata.SpacePointModuleLabel: @erase
+sbnd_gnewcalodata.NormTools: @local::sbnd_calonormtoolsdata  # YZ normalization for Data
 
 # gputnam 19Aug2024: Upate calibration to EMB
 # Values from: https://arxiv.org/abs/2407.12969