|
| 1 | + |
| 2 | +// Framework Includes |
| 3 | +#include "art/Framework/Core/EDProducer.h" |
| 4 | +#include "art/Framework/Principal/Event.h" |
| 5 | +#include "art/Framework/Principal/Handle.h" |
| 6 | +#include "art/Framework/Services/Registry/ServiceHandle.h" |
| 7 | +#include "art/Utilities/ToolMacros.h" |
| 8 | +#include "cetlib/cpu_timer.h" |
| 9 | +#include "fhiclcpp/ParameterSet.h" |
| 10 | +#include "messagefacility/MessageLogger/MessageLogger.h" |
| 11 | + |
| 12 | +#include "canvas/Persistency/Common/Ptr.h" |
| 13 | + |
| 14 | +// Tool include |
| 15 | +#include "larreco/Calorimetry/INormalizeCharge.h" |
| 16 | + |
| 17 | +// Services |
| 18 | +#include "lardata/DetectorInfoServices/DetectorClocksService.h" |
| 19 | + |
| 20 | +// ROOT includes |
| 21 | +#include "TFile.h" |
| 22 | +#include "TH2F.h" |
| 23 | + |
| 24 | +// C++ includes |
| 25 | +#include <map> |
| 26 | +#include <string> |
| 27 | +#include <stdexcept> |
| 28 | +#include <iostream> |
| 29 | + |
| 30 | + |
| 31 | + |
| 32 | +namespace sbnd { |
| 33 | + namespace calo { |
| 34 | + |
| 35 | + class NormalizeYZ : public INormalizeCharge{ |
| 36 | + public: |
| 37 | + explicit NormalizeYZ(fhicl::ParameterSet const& pset); |
| 38 | + |
| 39 | + void configure(const fhicl::ParameterSet& pset) override; |
| 40 | + |
| 41 | + double Normalize(double dQdx, |
| 42 | + const art::Event& evt, |
| 43 | + const recob::Hit& hit, |
| 44 | + const geo::Point_t& location, |
| 45 | + const geo::Vector_t& direction, |
| 46 | + double t0) override; |
| 47 | + |
| 48 | + private: |
| 49 | + void reconfigure(const fhicl::ParameterSet& pset); |
| 50 | + |
| 51 | + std::map<std::string, TH2F*> fCorrHists; |
| 52 | + std::string fFileName; |
| 53 | + bool fVerbose; |
| 54 | + }; |
| 55 | + |
| 56 | + NormalizeYZ::NormalizeYZ(fhicl::ParameterSet const& pset){ |
| 57 | + reconfigure(pset); // delegate config to reusable function |
| 58 | + } |
| 59 | + |
| 60 | + void NormalizeYZ::configure(const fhicl::ParameterSet& pset){ |
| 61 | + reconfigure(pset); |
| 62 | + } |
| 63 | + |
| 64 | + void NormalizeYZ::reconfigure(const fhicl::ParameterSet& pset){ |
| 65 | + fFileName = pset.get<std::string>("FileName"); |
| 66 | + fVerbose = pset.get<bool>("Verbose", false); |
| 67 | + |
| 68 | + std::string fname; |
| 69 | + cet::search_path sp("FW_SEARCH_PATH"); |
| 70 | + sp.find_file(fFileName, fname); |
| 71 | + |
| 72 | + TFile* f = TFile::Open(fname.c_str(), "READ"); |
| 73 | + if (!f || f->IsZombie()) { |
| 74 | + throw cet::exception("NormalizeYZ") << "Failed to open correction map file: " << fFileName; |
| 75 | + } |
| 76 | + |
| 77 | + for (int plane = 0; plane < 3; plane++){ // planes : 2 inductions (idx : 0, 1) and 1 collection (idx : 2) |
| 78 | + for (int tpc = 0; tpc < 2; tpc++){ // tpc : east (idx : 0) and west (idx : 1) TPCs |
| 79 | + std::string histname = Form("CzyHist_%d_%d", plane, tpc); |
| 80 | + TH2F* h = (TH2F*)f->Get(histname.c_str()); |
| 81 | + if (!h) { |
| 82 | + throw cet::exception("NormalizeYZ") << "Missing histogram in file: " << histname; |
| 83 | + } |
| 84 | + |
| 85 | + fCorrHists[Form("plane%d_%d", plane, tpc)] = (TH2F*)h->Clone(); |
| 86 | + fCorrHists[Form("plane%d_%d", plane, tpc)]->SetDirectory(nullptr); |
| 87 | + } |
| 88 | + } |
| 89 | + f->Close(); |
| 90 | + } |
| 91 | + |
| 92 | + double NormalizeYZ::Normalize(double dQdx, |
| 93 | + const art::Event& evt, |
| 94 | + const recob::Hit& hit, |
| 95 | + const geo::Point_t& location, |
| 96 | + const geo::Vector_t&, |
| 97 | + double){ |
| 98 | + |
| 99 | + int plane = hit.WireID().Plane; |
| 100 | + //int tpc = hit.WireID().TPC; |
| 101 | + |
| 102 | + // just to be sure and consistent with the logic of input YZ map |
| 103 | + // seems like current setup of directly calling hit.WireID().TPC has a tolerance of 30cm |
| 104 | + // meaning - an approx region of -30<x<30 lies in both the TPCs |
| 105 | + int tpc; |
| 106 | + if(location.X()<0){ |
| 107 | + tpc = 0; |
| 108 | + } else tpc = 1; |
| 109 | + |
| 110 | + |
| 111 | + std::string key = Form("plane%d_%d", plane, tpc); |
| 112 | + |
| 113 | + auto it = fCorrHists.find(key); |
| 114 | + if (it == fCorrHists.end()) { |
| 115 | + mf::LogWarning("NormalizeYZ") << "No correction histogram for " << key << ". Returning uncorrected dQdx"; |
| 116 | + return dQdx; |
| 117 | + } |
| 118 | + |
| 119 | + TH2F* hCorr = it->second; |
| 120 | + int binX = hCorr->GetXaxis()->FindBin(location.Z()); |
| 121 | + int binY = hCorr->GetYaxis()->FindBin(location.Y()); |
| 122 | + double scale = hCorr->GetBinContent(binX, binY); |
| 123 | + |
| 124 | + if (fVerbose){ |
| 125 | + std::cout << "[NormalizeYZ] Plane: " << plane << ", TPC: " << tpc |
| 126 | + << ", Y: " << location.Y() << ", Z: " << location.Z() |
| 127 | + << ", Scale: " << scale << ", dQdx (raw): " << dQdx |
| 128 | + << ", dQdx (corrected): " << dQdx / scale << std::endl; |
| 129 | + } |
| 130 | + |
| 131 | + if (scale < 1e-3){ |
| 132 | + mf::LogWarning("NormalizeYZ") << "Invalid scale at (Y,Z)=(" << location.Y() << "," << location.Z() << "). Returning uncorrected dQdx"; |
| 133 | + return dQdx; |
| 134 | + } |
| 135 | + |
| 136 | + return dQdx / scale; |
| 137 | + } |
| 138 | + |
| 139 | + DEFINE_ART_CLASS_TOOL(NormalizeYZ) |
| 140 | + |
| 141 | + } // namespace calo |
| 142 | +} // namespace sbnd |
| 143 | + |
0 commit comments