diff --git a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt index 96981ef6a..ccb73af6b 100644 --- a/sbncode/CAFMaker/RecoUtils/CMakeLists.txt +++ b/sbncode/CAFMaker/RecoUtils/CMakeLists.txt @@ -2,6 +2,7 @@ art_make_library( LIBRARY_NAME caf_RecoUtils SOURCE RecoUtils.cc LIBRARIES + sbnanaobj::StandardRecord art::Framework_Core art::Framework_Services_Registry art::Framework_Principal diff --git a/sbncode/HitFinder/BadChannelHitCheck_module.cc b/sbncode/HitFinder/BadChannelHitCheck_module.cc new file mode 100644 index 000000000..6ff63fe28 --- /dev/null +++ b/sbncode/HitFinder/BadChannelHitCheck_module.cc @@ -0,0 +1,127 @@ +//////////////////////////////////////////////////////////////////////// +// Class: BadChannelHitCheck +// Plugin Type: analyzer +// File: BadChannelHitCheck_module.cc +// +// Validation analyzer for the GaussHitFinderSBN "ExcludeBadChannels" +// feature. For each configured recob::Hit collection it checks every +// hit's channel against the channel-status database service; if any hit +// sits on a channel the database flags as bad, the art job is failed +// (by throwing). This lets a ctest verify end-to-end that bad channels +// are excluded from the hit finder output. +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/types/Atom.h" +#include "fhiclcpp/types/Sequence.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t +#include "lardataobj/RecoBase/Hit.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" + +#include + +namespace sbn { + + class BadChannelHitCheck : public art::EDAnalyzer { + public: + struct Config { + fhicl::Sequence HitLabels{ + fhicl::Name("HitLabels"), + fhicl::Comment("recob::Hit collections to validate against the bad-channel database")}; + fhicl::Atom RequirePresent{ + fhicl::Name("RequirePresent"), + fhicl::Comment("Also fail when a hit lands on a channel the database does not list as present"), + false}; + fhicl::Atom ThrowOnFailure{ + fhicl::Name("ThrowOnFailure"), + fhicl::Comment("Throw on the first bad-channel hit; if false, collect and throw in endJob()"), + true}; + }; + using Parameters = art::EDAnalyzer::Table; + + explicit BadChannelHitCheck(Parameters const& p); + + void analyze(art::Event const& e) override; + void endJob() override; + + private: + std::vector const fHitLabels; + bool const fRequirePresent; + bool const fThrowOnFailure; + + unsigned long fNHits{0}; ///< total hits checked + unsigned long fNBad{0}; ///< hits found on bad/not-present channels + bool fReportedStatus{false}; ///< whether the DB bad-channel count was logged + }; + + //------------------------------------------------------------------- + BadChannelHitCheck::BadChannelHitCheck(Parameters const& p) + : art::EDAnalyzer{p} + , fHitLabels{p().HitLabels()} + , fRequirePresent{p().RequirePresent()} + , fThrowOnFailure{p().ThrowOnFailure()} + {} + + //------------------------------------------------------------------- + void BadChannelHitCheck::analyze(art::Event const& e) + { + lariov::ChannelStatusProvider const& chanStatus = + art::ServiceHandle()->GetProvider(); + + if (!fReportedStatus) { + fReportedStatus = true; + mf::LogInfo("BadChannelHitCheck") + << "channel-status database reports " << chanStatus.BadChannels().size() + << " bad and " << chanStatus.NoisyChannels().size() << " noisy channel(s)."; + } + + for (art::InputTag const& tag : fHitLabels) { + auto const& hits = e.getProduct>(tag); + + for (recob::Hit const& hit : hits) { + raw::ChannelID_t const channel = hit.Channel(); + ++fNHits; + + bool const isBad = chanStatus.IsBad(channel); + bool const isAbsent = fRequirePresent && !chanStatus.IsPresent(channel); + if (!isBad && !isAbsent) continue; + + ++fNBad; + mf::LogError("BadChannelHitCheck") + << "recob::Hit on " << (isBad ? "BAD" : "NOT-PRESENT") << " channel " << channel + << " in product '" << tag.encode() << "', event " << e.id(); + + if (fThrowOnFailure) { + throw art::Exception(art::errors::LogicError) + << "BadChannelHitCheck: recob::Hit found on " + << (isBad ? "bad" : "not-present") << " channel " << channel + << " (product '" << tag.encode() << "', event " << e.id() << ")\n"; + } + } + } + } + + //------------------------------------------------------------------- + void BadChannelHitCheck::endJob() + { + mf::LogInfo("BadChannelHitCheck") + << "Checked " << fNHits << " hit(s); " << fNBad << " on bad/not-present channels."; + + if (fNBad > 0) { + throw art::Exception(art::errors::LogicError) + << "BadChannelHitCheck FAILED: " << fNBad << " of " << fNHits + << " hit(s) are on bad/not-present channels.\n"; + } + } + +} // namespace sbn + +DEFINE_ART_MODULE(sbn::BadChannelHitCheck) diff --git a/sbncode/HitFinder/CMakeLists.txt b/sbncode/HitFinder/CMakeLists.txt index 1d4870cef..dc1e13975 100644 --- a/sbncode/HitFinder/CMakeLists.txt +++ b/sbncode/HitFinder/CMakeLists.txt @@ -50,6 +50,7 @@ cet_build_plugin(GaussHitFinderSBN art::module ) cet_build_plugin(ChannelROIToWire art::module LIBRARIES ${MODULE_LIBRARIES}) cet_build_plugin(WireToChannelROI art::module LIBRARIES ${MODULE_LIBRARIES}) +cet_build_plugin(BadChannelHitCheck art::module LIBRARIES ${MODULE_LIBRARIES}) install_headers() diff --git a/sbncode/HitFinder/ChannelROIToWire_module.cc b/sbncode/HitFinder/ChannelROIToWire_module.cc index 42d304537..da7b1c4b4 100644 --- a/sbncode/HitFinder/ChannelROIToWire_module.cc +++ b/sbncode/HitFinder/ChannelROIToWire_module.cc @@ -141,7 +141,7 @@ void ChannelROIToWire::produce(art::Event& evt) std::vector dataVec(range.data().size()); - for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) dataVec[binIdx] = std::round(range.data()[binIdx] / ADCScaleFactor); + for(size_t binIdx = 0; binIdx < range.data().size(); binIdx++) dataVec[binIdx] = range.data()[binIdx] / ADCScaleFactor; ROIVec.add_range(startTick, std::move(dataVec)); } diff --git a/sbncode/HitFinder/GaussHitFinderSBN_module.cc b/sbncode/HitFinder/GaussHitFinderSBN_module.cc index c0a91cc91..8b404146d 100644 --- a/sbncode/HitFinder/GaussHitFinderSBN_module.cc +++ b/sbncode/HitFinder/GaussHitFinderSBN_module.cc @@ -45,6 +45,8 @@ #include "larcoreobj/SimpleTypesAndConstants/RawTypes.h" // raw::ChannelID_t #include "larreco/HitFinder/HitFilterAlg.h" #include "lardataobj/RecoBase/Hit.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusService.h" +#include "larevt/CalibrationDBI/Interface/ChannelStatusProvider.h" #include "sbnobj/ICARUS/TPC/ChannelROI.h" #include "sbncode/HitFinder/HitFinderUtilities/HitCreator.h" @@ -72,6 +74,7 @@ namespace hit { const bool fFilterHits; const bool fFillHists; + const bool fExcludeBadChannels; ///< skip input ChannelROIs on database-flagged bad channels const std::string fCalDataModuleLabel; const std::string fAllHitsInstanceName; @@ -99,6 +102,9 @@ namespace hit { //HitFilterAlg implementation is threadsafe. std::unique_ptr fHitFilterAlg; ///< algorithm used to filter out noise hits + /// Bad-channel database provider; non-null only when fExcludeBadChannels is true. + const lariov::ChannelStatusProvider* fChannelStatus = nullptr; + //only used when fFillHists is true and in single threaded mode. TH1F* fFirstChi2; TH1F* fChi2; @@ -111,6 +117,7 @@ namespace hit { : SharedProducer{pset} , fFilterHits(pset.get("FilterHits", false)) , fFillHists(pset.get("FillHists", false)) + , fExcludeBadChannels(pset.get("ExcludeBadChannels", false)) , fCalDataModuleLabel(pset.get("CalDataModuleLabel")) , fAllHitsInstanceName(pset.get("AllHitsInstanceName", "")) , fLongMaxHitsVec(pset.get>("LongMaxHits", std::vector() = {25, 25, 25})) @@ -127,6 +134,10 @@ namespace hit { , fPulseRatioCuts( pset.get>("PulseRatioCuts", std::vector() = {0.35, 0.40, 0.20})) { + // Cache the bad-channel database provider when channel exclusion is enabled. + if (fExcludeBadChannels) + fChannelStatus = &art::ServiceHandle()->GetProvider(); + if (fFillHists && art::Globals::instance()->nthreads() > 1u) { throw art::Exception(art::errors::Configuration) << "Cannot fill histograms when multiple threads configured, please set fFillHists to " @@ -223,6 +234,13 @@ namespace hit { auto const& wireReadoutGeom = art::ServiceHandle()->Get(); + // When bad-channel exclusion is enabled, snapshot the set of channels the + // status database flags as bad once here (single-threaded). The parallel + // loop below then cheaply skips any input ChannelROI on those channels. + const lariov::ChannelStatusProvider::ChannelSet_t badChannels = + fExcludeBadChannels ? fChannelStatus->BadChannels() + : lariov::ChannelStatusProvider::ChannelSet_t{}; + // ############################################### // ### Making a ptr vector to put on the event ### // ############################################### @@ -295,6 +313,9 @@ namespace hit { raw::ChannelID_t channel = channelROI->Channel(); + // Skip channels flagged as bad in the channel-status database. + if (fExcludeBadChannels && badChannels.count(channel)) return; + // get the WireID for this hit std::vector wids = wireReadoutGeom.ChannelToWire(channel); // for now, just take the first option returned from ChannelToWire diff --git a/sbncode/WireMod/Utility/WireModUtility.cc b/sbncode/WireMod/Utility/WireModUtility.cc index 126044753..87d5fefe6 100644 --- a/sbncode/WireMod/Utility/WireModUtility.cc +++ b/sbncode/WireMod/Utility/WireModUtility.cc @@ -108,6 +108,11 @@ std::vector> sys::WireModUtility::GetHitTa return target_roi_vec; } +//--- FillROIMatchedIDEMap --- +void sys::WireModUtility::FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset) +{ +} + //--- FillROIMatchedEdepMap --- void sys::WireModUtility::FillROIMatchedEdepMap(std::vector const& edepVec, std::vector const& wireVec, double offset) { @@ -249,11 +254,6 @@ std::mapE(); } - // calculate EDep properties by TrackID - std::map TrackIDMatchedPropertyMap; - for (auto const& track_edeps : TrackIDMatchedEDepMap) - TrackIDMatchedPropertyMap[track_edeps.first] = CalcPropertiesFromEdeps(track_edeps.second, offset); - // map it all out std::map> EDepMatchedSubROIMap; // keys are indexes of edepPtrVec, values are vectors of indexes of subROIPropVec std::map> TrackIDMatchedSubROIMap; // keys are TrackIDs, values are sets of indexes of subROIPropVec @@ -266,11 +266,18 @@ std::mapPositionToTPC(edep_ptr->MidPoint()); - const auto plane0 = wireReadout->FirstPlane(curTPCGeom.ID()); - double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - double zeroTick = detPropData.ConvertXToTicks(0, plane0.ID()); - auto edep_tick = ticksPercm * edep_ptr->X() + (zeroTick + offset) + tickOffset; - edep_tick = detPropData.ConvertXToTicks(edep_ptr->X(), plane0.ID()) + offset + tickOffset; + std::map planeProjEdepTick; + for (const auto& planeGeom : wireReadout->Iterate(curTPCGeom.ID())) + { + double projTick = detPropData.ConvertXToTicks(edep_ptr->X(), planeGeom.ID()) + offset + tickOffset; + readout::ROPID ropID = wireReadout->WirePlaneToROP(planeGeom.ID()); + if (planeProjEdepTick.count(ropID) == 1) + { + std::cout << "MatchEdepsToSubROIs - Warning, overwriting tick projection at" << ropID.toString() + << "\n Changing from " << planeProjEdepTick[ropID] << " to " << projTick << std::endl; + } + planeProjEdepTick[ropID] = projTick; + } // loop over subROIs unsigned int closest_hit = std::numeric_limits::max(); @@ -278,6 +285,8 @@ std::mapChannelToROP(subroi_prop.channel); + double edep_tick = planeProjEdepTick.at(subroiROP); if (edep_tick > subroi_prop.center-subroi_prop.sigma && edep_tick < subroi_prop.center+subroi_prop.sigma) { EDepMatchedSubROIMap[i_e].push_back(i_h); @@ -353,6 +362,12 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd { if (edep_ptr->StepLength() == 0) continue; + + geo::TPCID curTPCID; + try { + curTPCID = geometry->PositionToTPC(edep_ptr->MidPoint()).ID(); + } + catch(...) {continue;} // ignore non-active depositions edep_props.x = edep_ptr->X(); edep_props.y = edep_ptr->Y(); @@ -361,13 +376,11 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_props.dxdr = (edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_props.dydr = (edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_props.dzdr = (edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); - edep_props.dedr = edep_ptr->E() / edep_ptr->StepLength(); total_energy_all += edep_ptr->E(); - const geo::TPCGeo& curTPCGeom = geometry->PositionToTPC(edep_ptr->MidPoint()); - for (auto const& plane : wireReadout->Iterate(curTPCGeom.ID())) { + for (auto const& plane : wireReadout->Iterate(curTPCID)) { int i_p = plane.ID().Plane; auto scales = GetViewScaleValues(edep_props, plane.View()); scales_e_weighted[i_p].r_Q += edep_ptr->E()*scales.r_Q; @@ -409,7 +422,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = 0.; edep_col_properties.dydr = 0.; edep_col_properties.dzdr = 0.; - edep_col_properties.dedr = 0.; double total_energy = 0.0; @@ -428,7 +440,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr += edep_ptr->E()*(edep_ptr->EndX() - edep_ptr->StartX()) / edep_ptr->StepLength(); edep_col_properties.dydr += edep_ptr->E()*(edep_ptr->EndY() - edep_ptr->StartY()) / edep_ptr->StepLength(); edep_col_properties.dzdr += edep_ptr->E()*(edep_ptr->EndZ() - edep_ptr->StartZ()) / edep_ptr->StepLength(); - edep_col_properties.dedr += edep_ptr->E()*edep_ptr->E() / edep_ptr->StepLength(); } @@ -440,7 +451,6 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd edep_col_properties.dxdr = edep_col_properties.dxdr / total_energy; edep_col_properties.dydr = edep_col_properties.dydr / total_energy; edep_col_properties.dzdr = edep_col_properties.dzdr / total_energy; - edep_col_properties.dedr = edep_col_properties.dedr / total_energy; } @@ -451,18 +461,28 @@ sys::WireModUtility::TruthProperties_t sys::WireModUtility::CalcPropertiesFromEd } edep_col_properties.x_rms_noWeight = std::sqrt(edep_col_properties.x_rms_noWeight); - if (total_energy > 0) + if (total_energy > 0) { edep_col_properties.x_rms = std::sqrt(edep_col_properties.x_rms/total_energy); + } + + // get ticks, etc. if the deposition is active + geo::TPCID tpcGeomID; + try { + tpcGeomID = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}).ID(); + } + catch(...) {} - const geo::TPCGeo& tpcGeom = geometry->PositionToTPC({edep_col_properties.x, edep_col_properties.y, edep_col_properties.z}); - const auto plane0 = wireReadout->FirstPlane(tpcGeom.ID()); - double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now - edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; - edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; - edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; - edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; - edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; - edep_col_properties.total_energy = total_energy; + if (tpcGeomID.isValid) { + // projecting edep to plane0, so be mindful when accessing the tick value + const auto plane0 = wireReadout->FirstPlane(tpcGeomID); + double ticksPercm = detPropData.GetXTicksCoefficient(); // this should be by TPCID, but isn't building right now + edep_col_properties.tick = detPropData.ConvertXToTicks(edep_col_properties.x , plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_rms = ticksPercm*edep_col_properties.x_rms; + edep_col_properties.tick_rms_noWeight = ticksPercm*edep_col_properties.x_rms_noWeight; + edep_col_properties.tick_min = detPropData.ConvertXToTicks(edep_col_properties.x_min, plane0.ID()) + offset + tickOffset; + edep_col_properties.tick_max = detPropData.ConvertXToTicks(edep_col_properties.x_max, plane0.ID()) + offset + tickOffset; + edep_col_properties.total_energy = total_energy; + } return edep_col_properties; @@ -551,6 +571,20 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: if(temp_scale>0.001) scales.r_sigma *= temp_scale; } + if (applyXXWScale) + { + if (graph2Ds_Charge_XXW[plane] == nullptr || + graph2Ds_Sigma_XXW [plane] == nullptr ) + throw cet::exception("WireModUtility") + << "Tried to apply XXW scale factor, but could not find graphs. Check that you have set those in the utility."; + double thXW = ThetaXW(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ()); + temp_scale = graph2Ds_Charge_XXW[plane]->Interpolate(truth_props.x, thXW); + if(temp_scale>0.001) scales.r_Q *= temp_scale; + + temp_scale = graph2Ds_Sigma_XXW [plane]->Interpolate(truth_props.x, thXW); + if(temp_scale>0.001) scales.r_sigma *= temp_scale; + } + if (applyXZAngleScale) { if (splines_Charge_XZAngle[plane] == nullptr || @@ -570,7 +604,7 @@ sys::WireModUtility::ScaleValues_t sys::WireModUtility::GetViewScaleValues(sys:: scales.r_sigma *= splines_Sigma_YZAngle [plane]->Eval(ThetaYZ_PlaneRel(truth_props.dxdr, truth_props.dydr, truth_props.dzdr, plane_obj.ThetaZ())); } - if (applydEdXScale) + if(applydEdXScale) { if (splines_Charge_dEdX[plane] == nullptr || splines_Sigma_dEdX [plane] == nullptr ) @@ -635,6 +669,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, double q_orig = 0.0; double q_mod = 0.0; double scale_ratio = 1.0; + double sigma_distance = 0.0; // loop over the ticks for(size_t i_t = 0; i_t < roi_data.size(); ++i_t) @@ -643,6 +678,7 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, q_orig = 0.0; q_mod = 0.0; scale_ratio = 1.0; + sigma_distance = 0.0; // loop over the subs for (auto const& subroi_prop : subROIPropVec) @@ -652,12 +688,16 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, q_orig += gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); q_mod += gausFunc(i_t + roi_prop.begin, subroi_prop.center, scale_vals.r_sigma * subroi_prop.sigma, scale_vals.r_Q * subroi_prop.total_q); + sigma_distance += ((i_t + roi_prop.begin - subroi_prop.center)*(i_t + roi_prop.begin - subroi_prop.center) / (subroi_prop.sigma*subroi_prop.sigma))*\ + gausFunc(i_t + roi_prop.begin, subroi_prop.center, subroi_prop.sigma, subroi_prop.total_q); if (verbose) std::cout << " Incrementing q_orig by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << subroi_prop.sigma << ", " << subroi_prop.total_q << ")" << '\n' << " Incrementing q_mod by gausFunc(" << i_t + roi_prop.begin << ", " << subroi_prop.center << ", " << scale_vals.r_sigma * subroi_prop.sigma << ", " << scale_vals.r_Q * subroi_prop.total_q << ")" << std::endl; } + sigma_distance = sigma_distance / q_orig; + // do some sanity checks if (isnan(q_orig)) { @@ -668,6 +708,12 @@ void sys::WireModUtility::ModifyROI(std::vector & roi_data, if (verbose) std::cout << "WARNING: obtained q_mod = NaN... setting scale to 0" << std::endl; scale_ratio = 0.0; + } else if (q_orig < 0.01) { // check that this is a sane limit + if (verbose) std::cout << "WARNING: obtained q_orig < 0.01 ... setting scale to 1" << std::endl; + scale_ratio = 1.0; + } else if (sigma_distance > 9.) { + if (verbose) std::cout << "WARNING: sigma_distance > 3.. setting scale to 1" << std::endl; + scale_ratio = 1.0; } else { scale_ratio = q_mod / q_orig; } diff --git a/sbncode/WireMod/Utility/WireModUtility.hh b/sbncode/WireMod/Utility/WireModUtility.hh index 1ae74c152..4c9d9fb44 100644 --- a/sbncode/WireMod/Utility/WireModUtility.hh +++ b/sbncode/WireMod/Utility/WireModUtility.hh @@ -19,6 +19,7 @@ #include "lardataobj/RecoBase/Hit.h" #include "lardataobj/RecoBase/Wire.h" #include "lardataobj/Simulation/SimEnergyDeposit.h" +#include "lardataobj/Simulation/SimChannel.h" #include "larcoreobj/SimpleTypesAndConstants/PhysicalConstants.h" #include "nusimdata/SimulationBase/MCParticle.h" #include "nusimdata/SimulationBase/MCTruth.h" @@ -45,6 +46,7 @@ namespace sys { bool applyXXZAngleScale; // do we scale with X vs XZ angle? bool applyXdQdXScale; // do we scale with X vs dQ/dX? bool applyXZAngledQdXScale; // do we scale with XZ angle vs dQ/dX? + bool applyXXWScale; // do we scale with X vs ThXW? double readoutWindowTicks; // how many ticks are in the readout window? double tickOffset; // do we want an offset in the ticks? @@ -67,6 +69,8 @@ namespace sys { std::vector graph2Ds_Sigma_XdQdX; // the graphs for width correction in X vs dQ/dX std::vector graph2Ds_Charge_XZAngledQdX; // the graphs for charge correction in XZ angle vs dQ/dX std::vector graph2Ds_Sigma_XZAngledQdX; // the graphs for width correction in XZ angle vs dQ/dX + std::vector graph2Ds_Charge_XXW; // the graphs for the charge correction in XXW + std::vector graph2Ds_Sigma_XXW; // the graphs for the width correction in XXW // lets try making a constructor here // assume we can get a geometry service, a detector clcok, and a detector properties @@ -84,6 +88,7 @@ namespace sys { const bool& arg_ApplyXXZAngleScale = false, const bool& arg_ApplyXdQdXScale = false, const bool& arg_ApplyXZAngledQdXScale = false, + const bool& arg_ApplyXXWScale = true, const double& arg_TickOffset = 0) : geometry(geom), wireReadout(wireRead), @@ -97,13 +102,14 @@ namespace sys { applyXXZAngleScale(arg_ApplyXXZAngleScale), applyXdQdXScale(arg_ApplyXdQdXScale), applyXZAngledQdXScale(arg_ApplyXZAngledQdXScale), + applyXXWScale(arg_ApplyXXWScale), readoutWindowTicks(detProp.ReadOutWindowSize()), // the default A2795 (ICARUS TPC readout board) readout window is 4096 samples tickOffset(arg_TickOffset) // tick offset is for MC truth, default to zero and set only as necessary { } // typedefs - typedef std::pair ROI_Key_t; + typedef std::pair ROI_Key_t; typedef std::pair SubROI_Key_t; typedef struct ROIProperties @@ -145,8 +151,8 @@ namespace sys { float total_energy; float x_min; float x_max; - float tick_min; - float tick_max; + float tick_min; // On Plane0!! + float tick_max; // On Plane0!! float y; float z; double dxdr; @@ -158,8 +164,8 @@ namespace sys { } TruthProperties_t; // internal containers - std::map< ROI_Key_t,std::vector > ROIMatchedEdepMap; - std::map< ROI_Key_t,std::vector > ROIMatchedHitMap; + std::map< ROI_Key_t, std::vector > ROIMatchedEdepMap; + std::map< ROI_Key_t, std::vector > ROIMatchedHitMap; // some useful functions // geometries @@ -190,9 +196,8 @@ namespace sys { double ThetaXZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); //double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; // don't need to rotate Y for this angle double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; @@ -203,9 +208,8 @@ namespace sys { double ThetaYZ_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) { - double planeAngleRad = planeAngle * (util::pi() / 180.0); - double sinPlaneAngle = std::sin(planeAngleRad); - double cosPlaneAngle = std::cos(planeAngleRad); + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; @@ -213,8 +217,32 @@ namespace sys { double theta = std::atan2(dydrPlaneRel, dzdrPlaneRel); return FoldAngle(theta); } - // theste are set in the .cc file + + double ThetaXY_PlaneRel(double dxdr, double dydr, double dzdr, double planeAngle) + { + double sinPlaneAngle = std::sin(planeAngle); + double cosPlaneAngle = std::cos(planeAngle); + + double dydrPlaneRel = dydr * cosPlaneAngle - dzdr * sinPlaneAngle; + //double dzdrPlaneRel = dzdr * cosPlaneAngle + dydr * sinPlaneAngle; // don't need to rotate Z for this angle + + double theta = std::atan2(dxdr, dydrPlaneRel); + return FoldAngle(theta); + } + + double ThetaXW(double dxdr, double dydr, double dzdr, double planeAngle) + { + // Want angle to vert, so subtract pi + double sin = std::sin(planeAngle - ::util::pi<>()); + double cos = std::cos(planeAngle - ::util::pi<>()); + + double cosG = std::abs(dydr * sin + dzdr * cos); + double theta = std::atan(dxdr / cosG); + return std::abs(theta); + } + + // these are set in the .cc file ROIProperties_t CalcROIProperties(recob::Wire const&, size_t const&); std::vector> GetTargetROIs(sim::SimEnergyDeposit const&, double offset); @@ -222,6 +250,7 @@ namespace sys { void FillROIMatchedEdepMap(std::vector const&, std::vector const&, double offset); void FillROIMatchedHitMap(std::vector const&, std::vector const&); + void FillROIMatchedIDEMap(std::vector const& simchVec, std::vector const& wireVec, double offset); std::vector CalcSubROIProperties(ROIProperties_t const&, std::vector const&);