diff --git a/fcl/configurations/calibration_database_GlobalTags_icarus.fcl b/fcl/configurations/calibration_database_GlobalTags_icarus.fcl index 585af10c1..f67dd2296 100644 --- a/fcl/configurations/calibration_database_GlobalTags_icarus.fcl +++ b/fcl/configurations/calibration_database_GlobalTags_icarus.fcl @@ -5,7 +5,7 @@ BEGIN_PROLOG ICARUS_Calibration_GlobalTags: { - @table::TPC_CalibrationTags_Oct2023 + @table::TPC_CalibrationTags_Dec2024_TEMP_FOR_HARPS @table::PMT_CalibrationTags_Run2_Sept2023 @table::CRT_CalibrationTags_Oct2023 } diff --git a/fcl/configurations/calibration_database_TPC_TagSets_icarus.fcl b/fcl/configurations/calibration_database_TPC_TagSets_icarus.fcl index 056fdfe2f..e8715c38c 100644 --- a/fcl/configurations/calibration_database_TPC_TagSets_icarus.fcl +++ b/fcl/configurations/calibration_database_TPC_TagSets_icarus.fcl @@ -13,4 +13,14 @@ TPC_CalibrationTags_Oct2023: { } +## TEMPORARY +TPC_CalibrationTags_Dec2024_TEMP_FOR_HARPS: { + + tpc_channelstatus_data: "v3r1" + tpc_elifetime_data: "v2r1" + tpc_dqdxcalibration_data: "v2r1" + tpc_yz_correction_allplanes_data: "v1r0" + +} + END_PROLOG diff --git a/icaruscode/TPC/Calorimetry/NormalizeYZSQL_tool.cc b/icaruscode/TPC/Calorimetry/NormalizeYZSQL_tool.cc index 9c4226fcb..9868bd627 100644 --- a/icaruscode/TPC/Calorimetry/NormalizeYZSQL_tool.cc +++ b/icaruscode/TPC/Calorimetry/NormalizeYZSQL_tool.cc @@ -46,12 +46,14 @@ class NormalizeYZSQL : public INormalizeCharge class ScaleInfo { public: struct Point { + int iplane; int itpc; double y, z; }; class ScaleBin { public: + int iplane; int itpc; double ylo; double yhi; @@ -94,6 +96,7 @@ DEFINE_ART_CLASS_TOOL(NormalizeYZSQL) constexpr bool icarus::calo::NormalizeYZSQL::ScaleInfo::ScaleBin::operator< (const ScaleBin& other) const noexcept { + if (iplane != other.iplane) return iplane < other.iplane; if (itpc != other.itpc) return itpc < other.itpc; if (yhi != other.yhi) return yhi < other.yhi; return zhi < other.zhi; @@ -103,6 +106,7 @@ bool icarus::calo::NormalizeYZSQL::ScaleInfo::BinComp::operator() (const ScaleBin& a, const Point& b) const noexcept { // the bin `a` must be strictly before the point `b` + if (a.iplane != b.iplane) return a.iplane < b.iplane; if (a.itpc != b.itpc) return a.itpc < b.itpc; if (a.yhi <= b.y) return true; if (a.ylo > b.y) return false; @@ -113,6 +117,7 @@ bool icarus::calo::NormalizeYZSQL::ScaleInfo::BinComp::operator() (const Point& a, const ScaleBin& b) const noexcept { // the point `a` must be strictly before the bin `b` + if (a.iplane != b.iplane) return a.iplane < b.iplane; if (a.itpc != b.itpc) return a.itpc < b.itpc; if (a.y < b.ylo) return true; if (a.y >= b.yhi) return false; @@ -122,6 +127,7 @@ bool icarus::calo::NormalizeYZSQL::ScaleInfo::BinComp::operator() bool icarus::calo::NormalizeYZSQL::ScaleInfo::ScaleBin::contains (const Point& point) const noexcept { + if (point.iplane != iplane) return false; if (point.itpc != itpc) return false; if ((point.y < ylo) || (point.y >= yhi)) return false; return ((point.z >= zlo) && (point.z < zhi)); @@ -154,6 +160,7 @@ icarus::calo::NormalizeYZSQL::NormalizeYZSQL(fhicl::ParameterSet const &pset): void icarus::calo::NormalizeYZSQL::configure(const fhicl::ParameterSet& pset) {} const icarus::calo::NormalizeYZSQL::ScaleInfo& icarus::calo::NormalizeYZSQL::GetScaleInfo(uint64_t run) { + // check the cache if (fScaleInfos.count(run)) { return fScaleInfos.at(run); @@ -173,7 +180,13 @@ const icarus::calo::NormalizeYZSQL::ScaleInfo& icarus::calo::NormalizeYZSQL::Get // Iterate over the channels thisscale.bins.reserve(channels.size()); + for (unsigned ch = 0; ch < channels.size(); ch++) { + + long lplane; + fDB.GetNamedChannelData(ch, "plane", lplane); + int iplane(lplane); + std::string tpcname; fDB.GetNamedChannelData(ch, "tpc", tpcname); int itpc = -1; @@ -201,6 +214,7 @@ const icarus::calo::NormalizeYZSQL::ScaleInfo& icarus::calo::NormalizeYZSQL::Get bin.yhi = yhi; bin.zlo = zlo; bin.zhi = zhi; + bin.iplane = iplane; bin.itpc = itpc; bin.scale = scale; @@ -215,9 +229,13 @@ const icarus::calo::NormalizeYZSQL::ScaleInfo& icarus::calo::NormalizeYZSQL::Get double icarus::calo::NormalizeYZSQL::Normalize(double dQdx, const art::Event &e, const recob::Hit &hit, const geo::Point_t &location, const geo::Vector_t &direction, double t0) { + // Get the info ScaleInfo const& i = GetScaleInfo(e.id().runID().run()); + // plane + int plane = hit.WireID().Plane; + // compute itpc int cryo = hit.WireID().Cryostat; int tpc = hit.WireID().TPC; @@ -226,7 +244,7 @@ double icarus::calo::NormalizeYZSQL::Normalize(double dQdx, const art::Event &e, double y = location.y(); double z = location.z(); - ScaleInfo::Point const point { itpc, y, z }; + ScaleInfo::Point const point { plane, itpc, y, z }; ScaleInfo::ScaleBin const* b = i.findBin(point); double const scale = b? b->scale: 1; @@ -234,7 +252,7 @@ double icarus::calo::NormalizeYZSQL::Normalize(double dQdx, const art::Event &e, // TODO: what to do if no lifetime is found? throw an exception?? } - if (fVerbose) std::cout << "NormalizeYZSQL Tool -- Data Cryo: " << cryo << " TPC: " << tpc << " iTPC: " << itpc << " Y: " << y << " Z: " << z << " scale: " << scale << std::endl; + if (fVerbose) std::cout << "NormalizeYZSQL Tool -- Data Cryo: " << cryo << " Plane: " << plane << " TPC: " << tpc << " iTPC: " << itpc << " Y: " << y << " Z: " << z << " scale: " << scale << std::endl; return dQdx / scale; } diff --git a/icaruscode/TPC/Calorimetry/normtools_icarus.fcl b/icaruscode/TPC/Calorimetry/normtools_icarus.fcl index b873dd658..d755fd6eb 100644 --- a/icaruscode/TPC/Calorimetry/normtools_icarus.fcl +++ b/icaruscode/TPC/Calorimetry/normtools_icarus.fcl @@ -53,8 +53,8 @@ tpcgain_local: { yznorm_sql: { tool_type: NormalizeYZSQL - DBFileName: tpc_yz_correction_data - DBTag: @local::ICARUS_Calibration_GlobalTags.tpc_yz_correction_data + DBFileName: tpc_yz_correction_allplanes_data + DBTag: @local::ICARUS_Calibration_GlobalTags.tpc_yz_correction_allplanes_data Verbose: false } diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/CMakeLists.txt b/icaruscode/TPC/Tracking/ICARUSPandora/CMakeLists.txt index 48ab4da00..6bf8cd1a4 100644 --- a/icaruscode/TPC/Tracking/ICARUSPandora/CMakeLists.txt +++ b/icaruscode/TPC/Tracking/ICARUSPandora/CMakeLists.txt @@ -1,6 +1,43 @@ # where should the scripts/..xml file be installed? Perhaps in bin? +art_make_library( + LIBRARIES + lardataobj::RecoBase + fhiclcpp::fhiclcpp + cetlib::cetlib +) +set( MODULE_LIBRARIES + larcorealg::Geometry + larcorealg::Geometry + larcore::Geometry_Geometry_service + lardata::Utilities + larevt::Filters + lardataobj::RecoBase + lardata::ArtDataHelper + nurandom::RandomUtils_NuRandomService_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support ROOT::Core + art_root_io::TFileService_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + messagefacility::headers + fhiclcpp::fhiclcpp + cetlib::cetlib + CLHEP::Random + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + ROOT::FFTW + + ) +cet_build_plugin(HARPS art::module LIBRARIES ${MODULE_LIBRARIES}) + set( TOOL_LIBRARIES larpandora::LArPandoraInterface lardataobj::RecoBase @@ -12,6 +49,7 @@ set( TOOL_LIBRARIES cet_build_plugin(LArPandoraHitCollectionToolICARUS art::tool LIBRARIES ${TOOL_LIBRARIES}) +install_headers() install_fhicl() install_source() diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/HARPS_module.cc b/icaruscode/TPC/Tracking/ICARUSPandora/HARPS_module.cc new file mode 100644 index 000000000..93940991e --- /dev/null +++ b/icaruscode/TPC/Tracking/ICARUSPandora/HARPS_module.cc @@ -0,0 +1,638 @@ +//////////////////////////////////////////////////////////////////////// +// Class: HARPS (Hit Activity Removal from Particles for Systematics) +// Plugin Type: producer (Unknown Unknown) +// File: HARPS_module.cc +// +// Generated at Tue May 23 10:46:39 2023 by Bruce Howard using cetskelgen +// from version . +//////////////////////////////////////////////////////////////////////// + +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/Run.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/InputTag.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "canvas/Persistency/Common/FindMany.h" +#include "canvas/Persistency/Common/FindManyP.h" + +#include "cetlib_except/exception.h" + +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/RecoBase/Track.h" +#include "lardataobj/RecoBase/Vertex.h" +#include "lardataobj/RecoBase/Hit.h" +#include "lardata/ArtDataHelper/HitCreator.h" +#include "lardataobj/RecoBase/SpacePoint.h" + +// For random implementation: see +// https://github.com/SBNSoftware/icaruscode/blob/7b71f0213a1a66913cfe1cd795f49a9215570145/icaruscode/PMT/OpReco/FakePhotoS_module.cc#L16 +#include "nurandom/RandomUtils/NuRandomService.h" +#include "CLHEP/Random/RandFlat.h" + +#include +#include + +class HARPS; + + +class HARPS : public art::EDProducer { +public: + explicit HARPS(fhicl::ParameterSet const& p); + // The compiler-generated destructor is fine for non-base + // classes without bare pointers or other resource use. + + // Plugins should not be copied or assigned. + HARPS(HARPS const&) = delete; + HARPS(HARPS&&) = delete; + HARPS& operator=(HARPS const&) = delete; + HARPS& operator=(HARPS&&) = delete; + + // Required functions. + void produce(art::Event& e) override; + + double FindNearestTrajPointRR( const std::map, double>& trjPtToRR, const double* in_xyz ); + + typedef std::pair WireTimePair_t; + typedef std::pair PairOfWireTimePair_t; + +private: + + bool fOnlyOnePerEvent; + + std::vector fHitModuleLabels; + std::vector fPFParticleModuleLabels; + std::vector fTrackModuleLabels; + + std::string fInterestFileName; + + // TODO: how to handle the hits without spacepoints... + double fMaskBeginningDist; ///< How far from beginning of track to mask (in cm) + double fMaskEndDist; ///< How far from end of track to mask (in cm) + + // Allow for one to choose random groups of wires to skip on the different planes... + // TODO: Attempt made to handle things crossing boundaries, but possibly could be improved. + std::vector< unsigned int > fMaskNBreak; + std::vector< int > fMaskNWires; + + bool fTPCHitsWireAssn; + std::string fTPCHitCreatorInstanceName; + + std::map< std::string, std::vector > fParticleListCryo0; + std::map< std::string, std::vector > fParticleListCryo1; + + CLHEP::HepRandomEngine &fFlatEngine; + CLHEP::RandFlat *fFlatRand; ///< Random number generator as in PMT/OpReco/FakePhotoS_module.cc + + bool fTagDaughters; ///< Are we tagging daughter particles? + bool fKeepContext; ///< Are we doing the inverse operation (ie keeping only the context)? + bool fUseTrajPointForDist; ///< Are we making a map for track trajectory points to start/end distance to use +}; + + +HARPS::HARPS(fhicl::ParameterSet const& p) + : EDProducer{p}, + fOnlyOnePerEvent ( p.get< bool >("OnlyOnePerEvent", true) ), + fHitModuleLabels ( p.get< std::vector >("HitModuleLabels") ), + fPFParticleModuleLabels ( p.get< std::vector >("PFParticleModuleLabels") ), + fTrackModuleLabels ( p.get< std::vector >("TrackModuleLabels") ), + fInterestFileName ( p.get("InterestFileName") ), + fMaskBeginningDist ( p.get("MaskBeginningDist", -1.) ), + fMaskEndDist ( p.get("MaskEndDist", -1.) ), + fMaskNBreak ( p.get< std::vector >("MaskNBreak", {0}) ), + fMaskNWires ( p.get< std::vector >("MaskNWires", {-1}) ), + fTPCHitsWireAssn ( p.get< bool >("TPCHitsWireAssn", true) ), + fTPCHitCreatorInstanceName ( p.get("TPCHitCreatorInstanaceName","") ), + fFlatEngine ( art::ServiceHandle()->createEngine(*this, "HepJamesRandom", "Gen", p, "Seed") ), + fTagDaughters ( p.get< bool >("TagDaughters", false) ), + fKeepContext ( p.get< bool >("KeepContext", false) ), + fUseTrajPointForDist ( p.get< bool >("UseTrajPointForDist", false) ) +{ + if ( fPFParticleModuleLabels.size()!=fTrackModuleLabels.size() || fPFParticleModuleLabels.size()!=fHitModuleLabels.size() ) + throw cet::exception("HARPS") << "Error... InputTag vectors need to be same size..." << std::endl; // borrow from elsewhere + + if ( fMaskBeginningDist>0. && fMaskEndDist>0. ) + throw cet::exception("HARPS") << "Error... Should not try to mask the beginning and ALSO some distance from the end..." << std::endl; + + ////// + // TODO: there is something a little strange here with saving the daughters and keeping the end or start dist... FOR NOW, if + // we do fTagDaughters, do *NOT* allow fMaskBeginningDist or fMaskEndDist... and in this case the vertex we want to return + // is the parent particle's start/vertex position. + if ( (fKeepContext || fTagDaughters) && (fMaskBeginningDist>0. || fMaskEndDist>0. || fUseTrajPointForDist) ) + throw cet::exception("HARPS") << "Error... something strange about trying to keep daughters and remove based on start/end positions..." << std::endl; + + recob::HitCollectionCreator::declare_products(producesCollector(), fTPCHitCreatorInstanceName, fTPCHitsWireAssn, false); + + // Load in the PFP list of interest and put contents into a map for fast lookup. BASED ON larevt/Filters/EventFilter_module.cc + std::ifstream in; + in.open(fInterestFileName.c_str()); + char line[1024]; + while (1) { + in.getline(line, 1024); + if (!in.good()) break; + unsigned int nR, nS, nE, nC, nP; + sscanf(line, "%u %u %u %u %u", &nR, &nS, &nE, &nC, &nP); + std::string evtID = std::to_string(nR) + ":" + std::to_string(nS) + ":" + std::to_string(nE); + std::map< std::string, std::vector >& aliasParticleList = (nC == 0) ? fParticleListCryo0 : fParticleListCryo1; + if ( fOnlyOnePerEvent && (aliasParticleList.find(evtID) != aliasParticleList.end()) ) + continue; + aliasParticleList[evtID].push_back(nP); + } + in.close(); + + fFlatRand = new CLHEP::RandFlat(fFlatEngine,0.,1.); + + produces>(); +} + +void HARPS::produce(art::Event& e) +{ + // As in Gauss Hit Finder code (larreco/HitFinder/GausHitFinder_module.cc) + recob::HitCollectionCreator hitCol(e, fTPCHitCreatorInstanceName, fTPCHitsWireAssn, false); + + // Make the vector of vertices: if using fKeepContext then we will output an EMPTY vector of vertices + std::unique_ptr< std::vector > vtxVec( std::make_unique< std::vector >() ); + + // Return empty hit container if no particles of interest. + std::string evtID = std::to_string(e.run()) + ":" + std::to_string(e.subRun()) + ":" + std::to_string(e.event()); + if ( fParticleListCryo0.find(evtID) == fParticleListCryo0.end() && + fParticleListCryo1.find(evtID) == fParticleListCryo1.end() ) { + e.put(std::move(vtxVec)); + hitCol.put_into(e); + return; + } + + // DON'T THINK I ACTUALLY WANT THESE COMMENTS ANY MORE NOW THAT I FIGURED OUT THE ISSUE... + //std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + //std::cout << "~~HARPS~~ We're looking at event " << evtID << std::endl; + //std::cout << "~~HARPS~~ And we want PFP(s): East " << fParticleListCryo0[evtID].size() << ", West " << fParticleListCryo1[evtID].size() << std::endl; + //if ( fParticleListCryo0.find(evtID) != fParticleListCryo0.end() ) + // std::cout << "~~ Cryo0: " << fParticleListCryo0[evtID].size() << " ... " + // << (fParticleListCryo0[evtID].size() > 0 ? fParticleListCryo0[evtID][0] : 0) << std::endl; + //if ( fParticleListCryo1.find(evtID) == fParticleListCryo1.end() ) + // std::cout << "~~ Cryo1: " << fParticleListCryo1[evtID].size() << " ... " + // << (fParticleListCryo1[evtID].size() > 0 ? fParticleListCryo1[evtID][0] : 0) << std::endl; + //std::cout << "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!" << std::endl; + + //std::cout << "!!!! --- Got here 1" << std::endl; + + // Load in the PFParticles, Tracks, and Hits & the necessary FindManyP's + for ( unsigned int iCryo=0; iCryo<2; ++iCryo ) { + std::map< std::string, std::vector >& aliasParticleList = (iCryo == 0) ? fParticleListCryo0 : fParticleListCryo1; + if ( aliasParticleList.find(evtID) == aliasParticleList.end() ) continue; + + //std::cout << "Currently looking at Cryostat " << iCryo << " which has " << aliasParticleList[evtID].size() << " particles." << std::endl; + + art::InputTag prtLabel = (fPFParticleModuleLabels.size() > 1 ? fPFParticleModuleLabels[iCryo] : fPFParticleModuleLabels[0]); + art::InputTag trkLabel = (fTrackModuleLabels.size() > 1 ? fTrackModuleLabels[iCryo] : fTrackModuleLabels[0]); + art::InputTag hitLabel = (fHitModuleLabels.size() > 1 ? fHitModuleLabels[iCryo] : fHitModuleLabels[0]); + + // as in previous code + art::Handle< std::vector > pfpHandle; + std::vector< art::Ptr > pfps; + if ( e.getByLabel(prtLabel, pfpHandle) ) { + art::fill_ptr_vector(pfps, pfpHandle); + } + else { + mf::LogError("HARPS") << "Error pulling in PFParticle product... Skipping cryostat."; + } + + // How many/check + //std::cout << "!!!! --- Got here 2 --> there are " << pfps.size() << " pfps." << std::endl; + + // Get the find many for vertices, we only really want this if fTagDaughters is true at moment... + art::FindManyP fmVtx(pfpHandle, e, prtLabel); + if ( !fmVtx.isValid() ){ + throw cet::exception("HARPS") << "Module wants to use invalid association to recob::Vertex." << std::endl; + } + + // sort the PFParticles by id + // annoying as recob::PFParticle has a < opperator, but we have art::Ptr + std::sort(pfps.begin(), pfps.end(), [](art::Ptr pfpA, art::Ptr pfpB){ return pfpA->Self() < pfpB->Self(); }); + + // if we are tagging daughters find them here + if (fTagDaughters) + { + // if we're going to tag daughters, save the parent particles' vertices first + for ( auto const& pfpId : aliasParticleList[evtID] ) { + const std::vector< art::Ptr > vtxFromPFP = fmVtx.at(pfps[pfpId]->Self()); + if ( vtxFromPFP.size() != 1 ) continue; + double vtxPos[3] = {vtxFromPFP.at(0)->position().X(), vtxFromPFP.at(0)->position().Y(), vtxFromPFP.at(0)->position().Z()}; + vtxVec->push_back( recob::Vertex(vtxPos) ); + } + + auto idItr = aliasParticleList[evtID].begin(); + while (idItr != aliasParticleList[evtID].end()) + { + std::vector dghts = pfps[*idItr]->Daughters(); + aliasParticleList[evtID].insert(aliasParticleList[evtID].end(), dghts.begin(), dghts.end()); + ++idItr; + } + } + + // if we are keeping context, get the complement of our indices + if (fKeepContext) + { + // start with {0, 1, ..., pfps.size() - 1} + std::vector complement(pfps.size()); + std::iota(complement.begin(), complement.end(), 0); + + // make sure the initial pfp ids are sorted + // we want descending order so when we erase from complement we don't effect the indexing + std::sort(aliasParticleList[evtID].begin(), aliasParticleList[evtID].end(), [](size_t a, size_t b){ return a > b; }); + + // erase the indices in aliasParticleList[evtID] from complement + for (auto const& id : aliasParticleList[evtID]) + complement.erase(complement.begin() + id); + + // set the particle list to the complement + aliasParticleList[evtID] = complement; + } + + art::Handle< std::vector > trkHandle; + std::vector< art::Ptr > trks; + if ( e.getByLabel(trkLabel, trkHandle) ) { + art::fill_ptr_vector(trks, trkHandle); + } + else { + mf::LogError("HARPS") << "Error pulling in Track product... Skipping cryostat."; + } + + art::Handle< std::vector > hitHandle; + std::vector< art::Ptr > hits; + if ( e.getByLabel(hitLabel, hitHandle) ) { + art::fill_ptr_vector(hits, hitHandle); + } + else { + mf::LogError("HARPS") << "Error pulling in Hit product... Skipping cryostat."; + } + + // For the Wire Assns in HitCreator + art::FindManyP fmWire(hitHandle, e, hitLabel); + if ( !fmWire.isValid() && fTPCHitsWireAssn ){ + throw cet::exception("HARPS") << "Module wants to use Hit-Wire Associations, but fmWire invalid." << std::endl; + } + + art::FindManyP fmTrk(pfpHandle, e, trkLabel); + art::FindManyP fmHit(trkHandle, e, trkLabel); + art::FindManyP fmSps(hitHandle, e, prtLabel); + if ( !fmTrk.isValid() || !fmHit.isValid() || !fmSps.isValid() ){ + throw cet::exception("HARPS") << "Module wants to use invalid associations." << std::endl; + } + + //std::cout << "!!!! --- Got here 3 --> About to loop over the pfpIds... We have " + // << aliasParticleList[evtID].size() << " to care about." << std::endl; + + // Find the PFParticles we care about and put their hits into the collection + for (auto const& id : aliasParticleList[evtID]) { + art::Ptr pfp = pfps[id]; + + const std::vector< art::Ptr > trkFromPFP = fmTrk.at(pfp.key()); + if ( trkFromPFP.size() != 1 ) { + mf::LogError("HARPS") << "fmTrk gave us more (or less) than one track for this PFParticle... Skipping."; + continue; + } + + // If using Traj Points for distance calculation (via RR) + std::vector trkFilterTrajPt; + std::vector trkTrajPtLengths; + std::map, double> mapTrajPtEndDist; + + if ( fUseTrajPointForDist ) { + // First let's make a vector that is iterable (I think we can assume they are sorted?) + double prev_x(0.), prev_y(0.), prev_z(0.); + for ( unsigned int idxTrajPt=0; idxTrajPt < trkFromPFP.at(0)->NumberTrajectoryPoints(); ++idxTrajPt ) { + if ( !trkFromPFP.at(0)->HasValidPoint(idxTrajPt) ) continue; + trkFilterTrajPt.push_back( trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position ); + if ( idxTrajPt == 0 ) { + trkTrajPtLengths.push_back(0.); + prev_x = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.x(); + prev_y = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.y(); + prev_z = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.z(); + } + else { + double this_x = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.x(); + double this_y = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.y(); + double this_z = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position.z(); + trkTrajPtLengths.push_back( std::hypot( this_x-prev_x, this_y-prev_y, this_z-prev_z) ); + prev_x = this_x; + prev_y = this_y; + prev_z = this_z; + } + } + // Now let's build the maps + unsigned int idxTrajPt = 0; + for ( std::vector::iterator it=trkTrajPtLengths.begin(); it!=trkTrajPtLengths.end(); ++it ) { + const double rr = std::accumulate( it, trkTrajPtLengths.end(), 0. ); + //const recob::tracking::Point_t trajPtPoint = trkFromPFP.at(0)->TrajectoryPoint(idxTrajPt).position; + const std::tuple trajPtPos = {trkFilterTrajPt[idxTrajPt].x(),trkFilterTrajPt[idxTrajPt].y(),trkFilterTrajPt[idxTrajPt].z()}; + if ( mapTrajPtEndDist.find(trajPtPos) == mapTrajPtEndDist.end() ) { + mapTrajPtEndDist[trajPtPos] = rr; + } + else { + std::cout << "Found same traj point already? SKIPPING BUT PRINTING OUT THAT THIS HAS HAPPENED..." << std::endl; + } + idxTrajPt+=1; + } + } + + double trkStartX = trkFromPFP.at(0)->Start().X(); + double trkStartY = trkFromPFP.at(0)->Start().Y(); + double trkStartZ = trkFromPFP.at(0)->Start().Z(); + double trkEndX = trkFromPFP.at(0)->End().X(); + double trkEndY = trkFromPFP.at(0)->End().Y(); + double trkEndZ = trkFromPFP.at(0)->End().Z(); + + const std::vector< art::Ptr > hitsFromTrk = fmHit.at(trkFromPFP.at(0).key()); + + // Make a vector of the hits for each plane, and while doing it, compile a list of idx = Nw * 4 * 3 * Cryo + Nw * 3 * TPC + Nw * Plane + Wire, for now hardcode Nw = 5000 which should be long enough... + std::map< unsigned int, std::vector< art::Ptr > > planeToHitsMap; + std::map< unsigned int, std::set< unsigned int > > planeToIdxMap; + + // map for the interpolation + // first elem is PlaneID --> combination of cryo, tpc, plane + // second elem is pair of <, > + std::map< geo::PlaneID, PairOfWireTimePair_t > wireToTimeRangeMap; + + std::map< unsigned int, geo::PlaneID > planeNumToPlaneID_ForMaskBeginningDist; // planeID for begin point from which we'll search a box + std::map< unsigned int, WireTimePair_t > planeNumToWireTime_ForMaskBeginningDist; // WireTimePair_t for begin point from which we'll search a box + std::map< unsigned int, double > planeNumToPtDist_ForMaskBeginningDist; // min distance for begin point from which we'll search a box + std::map< unsigned int, geo::PlaneID > planeNumToPlaneID_ForMaskEndDist; // planeID for end point from which we'll search a box + std::map< unsigned int, WireTimePair_t > planeNumToWireTime_ForMaskEndDist; // WireTimePair_t for end point from which we'll search a box + std::map< unsigned int, double > planeNumToPtDist_ForMaskEndDist; // min distance for end point from which we'll search a box + + // Vertex coordinates: only used if fMaskBeginningDist or fMaskEndDist are > 0. + double minmaxDist = ( fMaskBeginningDist > 0. ? 9999. : ( fMaskEndDist > 0. ? -9999. : 0.) ); + double vtxX = -9999.; + double vtxY = -9999.; + double vtxZ = -9999.; + + // Loop as in my overlay module, copy hits passing the criteria into the maps for checking against other criteria... + for( auto const& iHitPtr : hitsFromTrk ) { + // Check if we want to remove this hit because we are masking X distance from the beginning: + // TODO: as a start, use the Pandora spacepoints<->hit association... is this safe or should we do something in 2d distances? + // UPDATE: loop through twice and save the max and min wires (and times?) the first time for things with a match and use these to + // save some of the 0 match cases the second time through via interpolation... only in the MaskBeginning or MaskEnd case... + if ( fMaskBeginningDist > 0. ) { + const std::vector< art::Ptr > spsFromHit = fmSps.at(iHitPtr.key()); + if ( spsFromHit.size() != 1 ) { + if ( spsFromHit.size() > 1 ) mf::LogError("HARPS") << "fmSps gave us something more than one spacepoint for this hit (" << spsFromHit.size() << ")... Skipping."; + continue; + } + const auto thisXYZ = spsFromHit.at(0)->XYZ(); + if ( std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ) < fMaskBeginningDist ) continue; + + // if this spacepoint is closest to the track start, save it as our 3D Vertex + if ( std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ) < minmaxDist ) { + minmaxDist = std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ); + vtxX = thisXYZ[0]; + vtxY = thisXYZ[1]; + vtxZ = thisXYZ[2]; + } + + // Check wireToTimeRangeMap + const auto planeid = iHitPtr->WireID().asPlaneID(); + const unsigned int thisplane = iHitPtr->WireID().Plane; + const unsigned int thiswire = iHitPtr->WireID().Wire; + const float thistime = iHitPtr->PeakTime(); + if ( wireToTimeRangeMap.count(planeid) == 0 ) { + wireToTimeRangeMap[planeid] = std::make_pair( std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()), + std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()) ); + } + else { + if ( thiswire < wireToTimeRangeMap[planeid].first.first ) { + wireToTimeRangeMap[planeid].first.first = thiswire; + wireToTimeRangeMap[planeid].first.second = thistime; + } + else if ( thiswire == wireToTimeRangeMap[planeid].first.first && thistime < wireToTimeRangeMap[planeid].first.second ) { + wireToTimeRangeMap[planeid].first.second = thistime; + } + else if ( thiswire == wireToTimeRangeMap[planeid].second.first && thistime > wireToTimeRangeMap[planeid].second.second ) { + wireToTimeRangeMap[planeid].second.second = thistime; + } + if ( thiswire > wireToTimeRangeMap[planeid].second.first ) { + wireToTimeRangeMap[planeid].second.first = thiswire; + wireToTimeRangeMap[planeid].second.second = thistime; + } + } // filling up the hit min/max map + + // Also deal with where to put a box to search near the beginning of track for hits unassociated with spacepoints + if ( planeNumToPtDist_ForMaskBeginningDist.count(thisplane) == 0 ) { + planeNumToPtDist_ForMaskBeginningDist[thisplane] = std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ); + planeNumToPlaneID_ForMaskBeginningDist[thisplane] = planeid; + planeNumToWireTime_ForMaskBeginningDist[thisplane] = std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()); + } + else if ( std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ) < planeNumToPtDist_ForMaskBeginningDist[thisplane] ) { + planeNumToPtDist_ForMaskBeginningDist[thisplane] = std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ); + planeNumToPlaneID_ForMaskBeginningDist[thisplane] = planeid; + planeNumToWireTime_ForMaskBeginningDist[thisplane] = std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()); + } + } // fMaskBeginningDist + else if ( fMaskEndDist > 0. ) { + const std::vector< art::Ptr > spsFromHit = fmSps.at(iHitPtr.key()); + if ( spsFromHit.size() != 1 ) { + if ( spsFromHit.size() > 1 ) mf::LogError("HARPS") << "fmSps gave us something more than one spacepoint for this hit (" << spsFromHit.size() << ")... Skipping."; + continue; + } + const auto thisXYZ = spsFromHit.at(0)->XYZ(); + + if ( !fUseTrajPointForDist ){ + // Just uses stright-line distance + if ( std::hypot( trkEndX-thisXYZ[0], trkEndY-thisXYZ[1], trkEndZ-thisXYZ[2] ) > fMaskEndDist ) continue; + + // if this spacepoint is furthest from the track end, save it as our 3D Vertex + if ( std::hypot( trkEndX-thisXYZ[0], trkEndY-thisXYZ[1], trkEndZ-thisXYZ[2] ) > minmaxDist ) { + minmaxDist = std::hypot( trkEndX-thisXYZ[0], trkEndY-thisXYZ[1], trkEndZ-thisXYZ[2] ); + vtxX = thisXYZ[0]; + vtxY = thisXYZ[1]; + vtxZ = thisXYZ[2]; + } + } + else { + // Uses the track's residual range by finding the closest traj point to this spacepoint + double nearestTrajPtRR = FindNearestTrajPointRR(mapTrajPtEndDist,thisXYZ); + if ( nearestTrajPtRR > fMaskEndDist || nearestTrajPtRR < 0. ) continue; + if ( nearestTrajPtRR > minmaxDist ) { + minmaxDist = nearestTrajPtRR; + vtxX = thisXYZ[0]; + vtxY = thisXYZ[1]; + vtxZ = thisXYZ[2]; + } + } + + // Check wireToTimeRangeMap + const auto planeid = iHitPtr->WireID().asPlaneID(); + const unsigned int thisplane = iHitPtr->WireID().Plane; + const unsigned int thiswire = iHitPtr->WireID().Wire; + const float thistime = iHitPtr->PeakTime(); + if ( wireToTimeRangeMap.count(planeid) == 0 ) { + wireToTimeRangeMap[planeid] = std::make_pair( std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()), + std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()) ); + } + else { + if ( thiswire < wireToTimeRangeMap[planeid].first.first ) { + wireToTimeRangeMap[planeid].first.first = thiswire; + wireToTimeRangeMap[planeid].first.second = thistime; + } + else if ( thiswire == wireToTimeRangeMap[planeid].first.first && thistime < wireToTimeRangeMap[planeid].first.second ) { + wireToTimeRangeMap[planeid].first.second = thistime; + } + else if ( thiswire == wireToTimeRangeMap[planeid].second.first && thistime > wireToTimeRangeMap[planeid].second.second ) { + wireToTimeRangeMap[planeid].second.second = thistime; + } + if ( thiswire > wireToTimeRangeMap[planeid].second.first ) { + wireToTimeRangeMap[planeid].second.first = thiswire; + wireToTimeRangeMap[planeid].second.second = thistime; + } + } // filling up the hit min/max map + + // Also deal with where to put a box to search near the beginning of track for hits unassociated with spacepoints + if ( planeNumToPtDist_ForMaskEndDist.count(thisplane) == 0 ) { + planeNumToPtDist_ForMaskEndDist[thisplane] = std::hypot( trkEndX-thisXYZ[0], trkEndY-thisXYZ[1], trkEndZ-thisXYZ[2] ); + planeNumToPlaneID_ForMaskEndDist[thisplane] = planeid; + planeNumToWireTime_ForMaskEndDist[thisplane] = std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()); + } + else if ( std::hypot( trkEndX-thisXYZ[0], trkEndY-thisXYZ[1], trkEndZ-thisXYZ[2] ) < planeNumToPtDist_ForMaskEndDist[thisplane] ) { + planeNumToPtDist_ForMaskEndDist[thisplane] = std::hypot( trkStartX-thisXYZ[0], trkStartY-thisXYZ[1], trkStartZ-thisXYZ[2] ); + planeNumToPlaneID_ForMaskEndDist[thisplane] = planeid; + planeNumToWireTime_ForMaskEndDist[thisplane] = std::make_pair(iHitPtr->WireID().Wire, iHitPtr->PeakTime()); + } + } // fMaskEndDist + + unsigned int plane = iHitPtr->WireID().Plane; + unsigned int idx = 5000*4*3*iHitPtr->WireID().Cryostat + 5000*3*iHitPtr->WireID().TPC + 5000*iHitPtr->WireID().Plane + iHitPtr->WireID().Wire; + planeToHitsMap[plane].push_back( iHitPtr ); + planeToIdxMap[plane].emplace( idx ); + } // loop hits + + // If we are masking by distance we may have skipped hits due to not being associated with a spacepoint... + // Let's try to add those back to the output by interpolating wires/times + if ( fMaskBeginningDist > 0. || fMaskEndDist > 0. ) { + for( auto const& iHitPtr : hitsFromTrk ) { + const std::vector< art::Ptr > spsFromHit = fmSps.at(iHitPtr.key()); + if ( spsFromHit.size() != 0 ) { + continue; // only want hits without spacepoints here + } + // The info to compare to map + const auto planeid = iHitPtr->WireID().asPlaneID(); + const unsigned int thisplane = iHitPtr->WireID().Plane; + const unsigned int thiswire = iHitPtr->WireID().Wire; + const float thistime = iHitPtr->PeakTime(); + if ( wireToTimeRangeMap.count(planeid) == 0 ) { + continue; // nothing on this plane in the map + } + // Check if it's within the interpolation and add it to the save map if so + if ( thiswire >= wireToTimeRangeMap[planeid].first.first && thistime >= wireToTimeRangeMap[planeid].first.second && + thiswire <= wireToTimeRangeMap[planeid].second.first && thistime <= wireToTimeRangeMap[planeid].second.second ) { + unsigned int plane = iHitPtr->WireID().Plane; + unsigned int idx = 5000*4*3*iHitPtr->WireID().Cryostat + 5000*3*iHitPtr->WireID().TPC + 5000*iHitPtr->WireID().Plane + iHitPtr->WireID().Wire; + planeToHitsMap[plane].push_back( iHitPtr ); + planeToIdxMap[plane].emplace( idx ); + } // check on interpolation + // Also check in a box near the beginning (if fMaskBeginningDist) or the end (if fMaskEndDist): let's start with ~3 cm + // --> so a box of about 10 wires and 50 ticks + else if ( ( fMaskBeginningDist > 0. && + planeid == planeNumToPlaneID_ForMaskBeginningDist[thisplane] && + abs(int(thiswire)-int(planeNumToWireTime_ForMaskBeginningDist[thisplane].first)) < 10 && + fabs(thistime-planeNumToWireTime_ForMaskBeginningDist[thisplane].second) < 50. ) || + ( fMaskEndDist > 0. && + planeid == planeNumToPlaneID_ForMaskEndDist[thisplane] && + abs(int(thiswire)-int(planeNumToWireTime_ForMaskEndDist[thisplane].first)) < 10 && + fabs(thistime-planeNumToWireTime_ForMaskEndDist[thisplane].second) < 50. ) ){ + unsigned int plane = iHitPtr->WireID().Plane; + unsigned int idx = 5000*4*3*iHitPtr->WireID().Cryostat + 5000*3*iHitPtr->WireID().TPC + 5000*iHitPtr->WireID().Plane + iHitPtr->WireID().Wire; + planeToHitsMap[plane].push_back( iHitPtr ); + planeToIdxMap[plane].emplace( idx ); + } // check "extrapolation box" + } // loop hits + + // PUT THE 3D VERTEX BACK INTO EVENT + double vtxPos[3] = {vtxX, vtxY, vtxZ}; + vtxVec->push_back( recob::Vertex(vtxPos) ); + } // only run if masking by distance + + // Pick wires to remove! + std::map< unsigned int, std::set< unsigned int > > planeToIdxRemoval; + if ( fMaskNWires.size() >= 1 && fMaskNBreak.size() >= 1 && fMaskNWires.size() == fMaskNBreak.size() ) { + + for ( unsigned int idxPlane = 0; idxPlane < fMaskNWires.size(); ++idxPlane ) { + if( fMaskNWires[idxPlane] < 0 || fMaskNBreak[idxPlane] == 0 ) continue; + if ( planeToIdxMap.find(idxPlane) == planeToIdxMap.end() ) continue; + std::vector planeIdxVec(planeToIdxMap[idxPlane].begin(),planeToIdxMap[idxPlane].end()); // vector copy for easy access... thanks Stack Overflow + unsigned int ithBreak = 0; + while ( ithBreak < fMaskNBreak[idxPlane] ) { + unsigned int locToBreak = planeToIdxMap[idxPlane].size()*fFlatRand->fire(0.,1.); + for ( unsigned int idxSubBreak=0; idxSubBreak < (unsigned int)fMaskNWires[idxPlane]; ++idxSubBreak ) { + if ( locToBreak+idxSubBreak >= planeToIdxMap[idxPlane].size() ) continue; + planeToIdxRemoval[idxPlane].emplace( planeIdxVec.at(locToBreak+idxSubBreak) ); + } // loop number of wires to break in this group + ithBreak+=1; + } // loop through number of groups of breaks to make in this plane + } // loop through number of planes on which to add breaks + } + + // Now loop over the hits we have and remove them if its wire didn't pass the test... + for ( auto const& [ plane, hitVec ] : planeToHitsMap ) { + for ( auto const& hit : hitVec ) { + unsigned int plane = hit->WireID().Plane; + unsigned int idx = 5000*4*3*hit->WireID().Cryostat + 5000*3*hit->WireID().TPC + 5000*hit->WireID().Plane + hit->WireID().Wire; + if ( planeToIdxRemoval.find(plane) != planeToIdxRemoval.end() && planeToIdxRemoval[plane].count( idx ) ) continue; + + // If we're not skipping this hit then write it to the file! + recob::Hit theHit = *hit; + if ( fTPCHitsWireAssn ) { + std::vector< art::Ptr > hitWires = fmWire.at(hit.key()); + if ( hitWires.size() == 0 ) throw cet::exception("HARPS") << "Hit found with no associated wires...\n"; + else if ( hitWires.size() > 1 ) mf::LogWarning("HARPS") << "Hit with >1 recob::Wire associated..."; + hitCol.emplace_back(theHit,hitWires[0]); + } + else hitCol.emplace_back(theHit); + } // loop hits + } // loop planes with hits + } // loop pfps + } // loop cryos + + e.put(std::move(vtxVec)); + hitCol.put_into(e); +} + +double HARPS::FindNearestTrajPointRR( const std::map, double>& trjPtToRR, const double* in_xyz ) +{ + //std::cout << "The hit spacepoint at (" << in_xyz[0] << ", " << in_xyz[1] << ", " << in_xyz[2] << ") "; + + if ( trjPtToRR.empty() ) { + // std::cout << " ... has no match." << std::endl; + return -5.; + } + + double minDistToPt = std::numeric_limits::max(); + double rrAtMinDist = -5.; + //double xAtMinDist = -5.; + //double yAtMinDist = -5.; + //double zAtMinDist = -5.; + + for ( auto const& [pt, rr] : trjPtToRR ) { + double ptx = std::get<0>(pt); + double pty = std::get<1>(pt); + double ptz = std::get<2>(pt); + if ( std::hypot(in_xyz[0]-ptx, in_xyz[1]-pty, in_xyz[2]-ptz) < minDistToPt ) { + minDistToPt = std::hypot(in_xyz[0]-ptx, in_xyz[1]-pty, in_xyz[2]-ptz); + rrAtMinDist = rr; + //xAtMinDist = ptx; + //yAtMinDist = pty; + //zAtMinDist = ptz; + } + } + + //std::cout << " ... has matched traj point at (" << xAtMinDist << ", " << yAtMinDist << ", " << zAtMinDist << "), with rr " << rrAtMinDist << std::endl; + + return rrAtMinDist; +} + + +DEFINE_ART_MODULE(HARPS) diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Cosmic_ICARUS_HARPS.xml b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Cosmic_ICARUS_HARPS.xml new file mode 100644 index 000000000..76fa5a9c4 --- /dev/null +++ b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Cosmic_ICARUS_HARPS.xml @@ -0,0 +1,324 @@ + + + false + true + true + + + + CaloHitListU + CaloHitListV + CaloHitListW + CaloHitList2D + CaloHitList2D + + + + + + CaloHitListU + ClustersU + true + true + + + + + + + + + + true + true + + + + + + + + + + + CaloHitListV + ClustersV + true + true + + + + + + + + + + true + true + + + + + + + + + + + CaloHitListW + ClustersW + true + true + + + + + + + + + + true + true + + + + + + + + + + + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + true + true + false + false + + + 0.750.75 + 0.750.75 + + + + + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + + + + 10. + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + + + + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + + false + ClustersU ClustersV ClustersW + + + + + + CaloHitListW + ClustersW + false + false + + + + CaloHitListU + ClustersU + false + false + + + + CaloHitListV + ClustersV + false + false + + + + + MuonParticles3D + DeltaRayParticles3D + + + ClustersW + MuonParticles3D + DeltaRayParticles3D + + + ClustersU + MuonParticles3D + DeltaRayParticles3D + + + ClustersV + MuonParticles3D + DeltaRayParticles3D + + + MuonParticles3D + DeltaRayParticles3D + ClustersU + ClustersV + ClustersW + + + DeltaRayParticles3D + + + + + ClustersW + + + ClustersU + + + ClustersV + + + ClustersW + + + ClustersU + + + ClustersV + + + + + 5 + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + + + + 2 + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + + + + + + 1 + ClustersU + ClustersV + ClustersW + MuonParticles3D + + + MuonParticles3D DeltaRayParticles3D + ClustersU ClustersV ClustersW + false + + + + + MuonParticles3D + MuonCaloHits3D + MuonClusters3D + + 3 + 3 + 3 + 3 + 2 + 2 + 2 + + + + MuonParticles3D + MuonCaloHits3D + MuonClusters3D + + + 5. + 3. + + + 5. + 2. + + + + + DeltaRayParticles3D + DeltaRayCaloHits3D + DeltaRayClusters3D + + 5.3. + 5.2. + + + + + + + MuonParticles3D + CRVertices3D + + + + + MuonParticles3D + CRVertices3D + + + + + MuonParticles3D DeltaRayParticles3D + CRVertices3D + ClustersU ClustersV ClustersW MuonClusters3D DeltaRayClusters3D + CaloHitListU CaloHitListV CaloHitListW CaloHitList2D + MuonParticles3D + + diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Master_ICARUS_HARPS.xml b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Master_ICARUS_HARPS.xml new file mode 100644 index 000000000..2adda634f --- /dev/null +++ b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Master_ICARUS_HARPS.xml @@ -0,0 +1,62 @@ + + + false + true + true + + + 1000. + + + + + CaloHitListU + CaloHitListV + CaloHitListW + CaloHitListCustom + CaloHitList2D + CaloHitList2D + + + CaloHitListU CaloHitListV CaloHitListW + true + + + + PandoraSettings_Cosmic_ICARUS_HARPS.xml + PandoraSettings_Neutrino_ICARUS_HARPS.xml + PandoraSettings_Slicing_ICARUS_HARPS.xml + + true + false + + + + 2.0 + cautious + + + + + PandoraBdt_v08_33_00_SBND.xml + NeutrinoId + 0 + 999 + true + + + Input + Input + false + RecreatedPfos + RecreatedClusters + RecreatedVertices + false + 2.0 + + + + true + true + + diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Neutrino_ICARUS_HARPS.xml b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Neutrino_ICARUS_HARPS.xml new file mode 100644 index 000000000..e24526e2b --- /dev/null +++ b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Neutrino_ICARUS_HARPS.xml @@ -0,0 +1,350 @@ + + + + true + true + true + + + + CaloHitListU + CaloHitListV + CaloHitListW + CaloHitListCustom + CaloHitList2D + CaloHitList2D + + + + + + CaloHitListU + ClustersU + true + true + + + + + + + + + + + + + + + + + + CaloHitListV + ClustersV + true + true + + + + + + + + + + + + + + + + + + CaloHitListW + ClustersW + true + true + + + + + + + + + + + + + + + + + + ClustersU ClustersV ClustersW + + + CaloHitListCustom + NeutrinoVertices3D + true + + + ClustersU ClustersV ClustersW + true + + + ClustersU ClustersV ClustersW + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + true + true + false + false + + + 0.750.75 + 0.750.75 + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + 5. + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + + + TrackParticles3D + ShowerParticles3D + false + + + ShowerParticles3D + + + ClustersU ClustersV ClustersW + true + + + ClustersU ClustersV ClustersW + + + ClustersU + ClustersV + ClustersW + ShowerParticles3D + + + + + + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + true + true + false + false + + + 0.750.75 + 0.750.75 + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + 5. + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + + + ClustersU ClustersV ClustersW + TrackParticles3D + + + ClustersU ClustersV ClustersW + TrackParticles3D + + + ClustersU ClustersV ClustersW + TrackParticles3D + true + 0.5 + 5 + 1. + + + + + ShowerParticles3D + ClustersU ClustersV ClustersW + + + ShowerParticles3D + ClustersU ClustersV ClustersW + + + ShowerParticles3D + ClustersU ClustersV ClustersW + + + + + TrackParticles3D + ShowerParticles3D + true + false + + + TrackParticles3D + TrackCaloHits3D + TrackClusters3D + + 3 + 3 + 3 + 3 + 2 + 2 + 2 + + + + ShowerParticles3D + ShowerCaloHits3D + ShowerClusters3D + + + + + + + + + + TrackParticles3D ShowerParticles3D + ClustersU ClustersV ClustersW TrackClusters3D ShowerClusters3D + + + TrackParticles3D ShowerParticles3D + ClustersU ClustersV ClustersW + + + TrackParticles3D ShowerParticles3D + ClustersU ClustersV ClustersW + true + + + + + NeutrinoVertices3D + NeutrinoParticles3D + + + NeutrinoParticles3D + TrackParticles3D ShowerParticles3D + false + + + + + + + + NeutrinoParticles3D + DaughterVertices3D + + + TrackParticles3D + ShowerParticles3D + true + true + PandoraBdt_v08_33_00_SBND.xml + PfoCharBDT + PandoraBdt_v08_33_00_SBND.xml + PfoCharBDTNoChargeInfo + + + + + + + + + + + + + + + + NeutrinoParticles3D + + + + + TrackParticles3D + DaughterVertices3D + + + + + NeutrinoParticles3D TrackParticles3D ShowerParticles3D + NeutrinoVertices3D DaughterVertices3D CandidateVertices3D + ClustersU ClustersV ClustersW TrackClusters3D ShowerClusters3D + CaloHitListU CaloHitListV CaloHitListW CaloHitList2D + NeutrinoParticles3D + + diff --git a/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Slicing_ICARUS_HARPS.xml b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Slicing_ICARUS_HARPS.xml new file mode 100644 index 000000000..830449897 --- /dev/null +++ b/icaruscode/TPC/Tracking/ICARUSPandora/scripts/PandoraSettings_Slicing_ICARUS_HARPS.xml @@ -0,0 +1,266 @@ + + + true + false + true + + + + CaloHitListU + CaloHitListV + CaloHitListW + CaloHitListCustom + CaloHitList2D + CaloHitList2D + + + + + + CaloHitListU + ClustersU + true + true + + + + + + + + + + + + + + + + + + CaloHitListV + ClustersV + true + true + + + + + + + + + + + + + + + + + + CaloHitListW + ClustersW + true + true + + + + + + + + + + + + + + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + true + true + false + false + + + 0.750.75 + 0.750.75 + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + 5. + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + + + TrackParticles3D + ShowerParticles3D + false + + + ShowerParticles3D + + + ClustersU ClustersV ClustersW + true + + + ClustersU ClustersV ClustersW + + + ClustersU + ClustersV + ClustersW + ShowerParticles3D + + + + + + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + true + true + false + false + + + 0.750.75 + 0.750.75 + + + + + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + 5. + ClustersU + ClustersV + ClustersW + TrackParticles3D + + + + + + + + + ClustersU ClustersV ClustersW + TrackParticles3D + + + ClustersU ClustersV ClustersW + TrackParticles3D + + + ClustersU ClustersV ClustersW + TrackParticles3D + true + 0.5 + 5 + 1. + + + + + TrackParticles3D + ShowerParticles3D + false + true + + + TrackParticles3D + TrackCaloHits3D + TrackClusters3D + + 3 + 3 + 3 + 3 + 2 + 2 + 2 + + + + ShowerParticles3D + ShowerCaloHits3D + ShowerClusters3D + + + + + + + + + + CaloHitListU + CaloHitListV + CaloHitListW + CaloHitListCustom + ClustersU + ClustersV + ClustersW + SliceClusters + SliceParticles + + TrackParticles3D + ShowerParticles3D + + + TrackParticles3D ShowerParticles3D + ClustersU ClustersV ClustersW TrackClusters3D ShowerClusters3D + + + + SliceParticles + +