From 2807b55d45a732804d4d868ecf9510dc2ab1b729 Mon Sep 17 00:00:00 2001 From: Felipe Wieler Date: Tue, 18 Nov 2025 08:39:38 -0600 Subject: [PATCH 1/2] Adding CVN code --- icaruscode/ICARUSCVN/CMakeLists.txt | 6 + icaruscode/ICARUSCVN/README.md | 77 ++ icaruscode/ICARUSCVN/fcls/CMakeLists.txt | 2 + icaruscode/ICARUSCVN/fcls/ICARUSCVNMapper.fcl | 188 ++++ icaruscode/ICARUSCVN/fcls/getTrueInfo.py | 74 ++ .../ICARUSCVN/fcls/run_ICARUSCVNEvaluator.fcl | 62 ++ .../ICARUSCVN/fcls/run_ICARUSCVNMapper.fcl | 72 ++ .../ICARUSCVN/fcls/run_ICARUSZlibMaker.fcl | 45 + .../ICARUSCVN/module_helpers/CMakeLists.txt | 43 + .../module_helpers/ICARUSICVNMapper.cxx | 142 +++ .../module_helpers/ICARUSICVNMapper.h | 61 ++ .../module_helpers/ICARUSICVNZlibMaker.cxx | 52 + .../module_helpers/ICARUSICVNZlibMaker.h | 87 ++ .../module_helpers/ICARUSITFNetHandler.h | 29 + .../module_helpers/ICARUSPixelMapProducer.cxx | 272 +++++ .../module_helpers/ICARUSPixelMapProducer.h | 49 + icaruscode/ICARUSCVN/module_helpers/classes.h | 8 + .../ICARUSCVN/module_helpers/classes_def.xml | 23 + icaruscode/ICARUSCVN/modules/CMakeLists.txt | 100 ++ .../modules/ICARUSCVNEvaluator_module.cc | 338 ++++++ .../modules/ICARUSCVNMapper_module.cc | 17 + .../modules/ICARUSCVNZlibMaker_module.cc | 999 ++++++++++++++++++ .../modules/ICARUSTFNetHandler_tool.cc | 298 ++++++ icaruscode/ICARUSCVN/tf/CMakeLists.txt | 9 + icaruscode/ICARUSCVN/tf/tf_graph.cc | 423 ++++++++ icaruscode/ICARUSCVN/tf/tf_graph.h | 80 ++ 26 files changed, 3556 insertions(+) create mode 100644 icaruscode/ICARUSCVN/CMakeLists.txt create mode 100644 icaruscode/ICARUSCVN/README.md create mode 100644 icaruscode/ICARUSCVN/fcls/CMakeLists.txt create mode 100644 icaruscode/ICARUSCVN/fcls/ICARUSCVNMapper.fcl create mode 100644 icaruscode/ICARUSCVN/fcls/getTrueInfo.py create mode 100644 icaruscode/ICARUSCVN/fcls/run_ICARUSCVNEvaluator.fcl create mode 100644 icaruscode/ICARUSCVN/fcls/run_ICARUSCVNMapper.fcl create mode 100644 icaruscode/ICARUSCVN/fcls/run_ICARUSZlibMaker.fcl create mode 100644 icaruscode/ICARUSCVN/module_helpers/CMakeLists.txt create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.cxx create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.h create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.cxx create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.h create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.cxx create mode 100644 icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h create mode 100644 icaruscode/ICARUSCVN/module_helpers/classes.h create mode 100644 icaruscode/ICARUSCVN/module_helpers/classes_def.xml create mode 100644 icaruscode/ICARUSCVN/modules/CMakeLists.txt create mode 100644 icaruscode/ICARUSCVN/modules/ICARUSCVNEvaluator_module.cc create mode 100644 icaruscode/ICARUSCVN/modules/ICARUSCVNMapper_module.cc create mode 100644 icaruscode/ICARUSCVN/modules/ICARUSCVNZlibMaker_module.cc create mode 100644 icaruscode/ICARUSCVN/modules/ICARUSTFNetHandler_tool.cc create mode 100644 icaruscode/ICARUSCVN/tf/CMakeLists.txt create mode 100644 icaruscode/ICARUSCVN/tf/tf_graph.cc create mode 100644 icaruscode/ICARUSCVN/tf/tf_graph.h diff --git a/icaruscode/ICARUSCVN/CMakeLists.txt b/icaruscode/ICARUSCVN/CMakeLists.txt new file mode 100644 index 000000000..73c8ebfcc --- /dev/null +++ b/icaruscode/ICARUSCVN/CMakeLists.txt @@ -0,0 +1,6 @@ +add_subdirectory(module_helpers) +add_subdirectory(tf) +add_subdirectory(modules) +add_subdirectory(fcls) +install_fhicl() + diff --git a/icaruscode/ICARUSCVN/README.md b/icaruscode/ICARUSCVN/README.md new file mode 100644 index 000000000..d63efcb73 --- /dev/null +++ b/icaruscode/ICARUSCVN/README.md @@ -0,0 +1,77 @@ +ICARUS CVN Mapper and Evaluator +Author: Felipe A. Wieler +Date: November 12, 2025 +Email: f3aw@proton.me + +------------------------------------------------------------ +Overview +------------------------------------------------------------ +The ICARUSCVNMapper and ICARUSCVNEvaluator modules provide the interface between reconstructed slices and the Convolutional Visual Network (CVN) models used in ICARUS. + +- ICARUSCVNMapper.fcl: Configuration file for detector mapping and label management, including Grid Configuration. +- ICARUSZlibMaker.fcl: Converts the mapper output into compressed pixel map archives. +- ICARUSCVNEvaluator.fcl: Runs the trained CVN TensorFlow models to evaluate each slice and produce classification scores. + +These tools together enable pixel map creation, compression, and evaluation for neutrino flavors classification in ICARUS data. + +------------------------------------------------------------ +Prerequisites +------------------------------------------------------------ +Before running, ensure the following: +- LArSoft and ICARUS code are correctly set up in your environment. +- TensorFlow model files are accessible, those are defined on the ICARUSCVNMapper.fcl. +- Input files (stage1) have been produced by standard ICARUS reconstruction chains. +- standard_cvnzlibmaker is properly configured, including the OutputDir path also in ICARUSCVNMapper.fcl. + +------------------------------------------------------------ +Pixel Map Generation +------------------------------------------------------------ +To create pixel maps from a stage1 file: +lar -c run_ICARUSCVNMapper.fcl -s + +This produces a ROOT output with slice and pixel map information. + +Next, run the zlib maker to generate the compressed maps: +lar -c run_ICARUSZlibMaker.fcl -s + +This step creates: +- .gz files → contain pixel map arrays +- .info files → include truth and slice metadata from Pandora + +These files are stored under the directory specified by OutputDir in standard_cvnzlibmaker on ICARUSCVNMapper.fcl. + +------------------------------------------------------------ +Running the Evaluator +------------------------------------------------------------ +To evaluate slices using the trained CVN models: +lar -c run_ICARUSCVNEvaluator.fcl -s + +This produces a new ROOT file containing the CryoE and CryoW classification scores for each slice. + +------------------------------------------------------------ +Output Summary +------------------------------------------------------------ +File Type | Description +----------|------------------------------------------------ +.gz | Compressed pixel map arrays +.info | Truth and metadata from Pandora per slice +.root | Evaluation results including CryoE and CryoW scores + +------------------------------------------------------------ +Example Workflow +------------------------------------------------------------ + +Pixel Map Creation: +lar -c run_ICARUSCVNMapper.fcl -s reco_stage1.root +lar -c run_ICARUSZlibMaker.fcl -s output_mapper.root +------------------------------------------------------------ + +Root Evaluation: +lar -c run_ICARUSCVNEvaluator.fcl -s reco_stage1.root + +------------------------------------------------------------ +Notes +------------------------------------------------------------ +- Use consistent stage1 files across all steps to maintain slice–truth consistency. +- Evaluation requires access to trained CVN TensorFlow models. +- The modules support both Cryostat East (CryoE) and Cryostat West (CryoW) toghether or individualy by changin the corresponding FHiCL file. diff --git a/icaruscode/ICARUSCVN/fcls/CMakeLists.txt b/icaruscode/ICARUSCVN/fcls/CMakeLists.txt new file mode 100644 index 000000000..fb3c8a297 --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/CMakeLists.txt @@ -0,0 +1,2 @@ +install_fhicl() + diff --git a/icaruscode/ICARUSCVN/fcls/ICARUSCVNMapper.fcl b/icaruscode/ICARUSCVN/fcls/ICARUSCVNMapper.fcl new file mode 100644 index 000000000..c9165eba8 --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/ICARUSCVNMapper.fcl @@ -0,0 +1,188 @@ +BEGIN_PROLOG + +# This file defines the configuration for the CVNMapper and CVNMapperProtoDUNE modules +standard_pixelmapproducer: +{ + TdcWidth: 640 #400 # 1000 + WireLength: 6500 #3064 #2205 #6223 #8500 #or 9000 #collection plane : 2205, induction plane : 6223 + TimeResolution: 1500 #1600 # 2560 # 3400 + Threshold: 0 + MultipleDrifts: true + verbose: false + ChangeWireNo: true + ReadoutSize: 4096. # in time ticks # 2560 + ShiftT: 850 #850. # 420. #size of the back/front porch + InductionWires: 6500 #3064 #6223 #2205 #1984 #number of wires in the first induction plane + FlipInductionView: false + UseT: true #use time 0 of pfp +} + +standard_cvnmapper_west: +{ + module_type: ICARUSCVNMapper + HitsModuleLabel: "cluster3DCryoW" + ClusterPMLabel: "cvnmap" + MinClusterHits: 100 + verbose: true + UseSlice: true + SliceLabel: "pandoraGausCryoW" + PFParticleModuleLabel: "pandoraGausCryoW" + T0Label: "pandoraGausCryoW" + PixelMapProducer: @local::standard_pixelmapproducer + MapVecSize: 100000 + +} + +standard_cvnmapper_east: +{ + module_type: ICARUSCVNMapper + HitsModuleLabel: "cluster3DCryoE" + ClusterPMLabel: "cvnmap" + MinClusterHits: 100 + verbose: true + UseSlice: true + SliceLabel: "pandoraGausCryoE" + PFParticleModuleLabel: "pandoraGausCryoE" + T0Label: "pandoraGausCryoE" + PixelMapProducer: @local::standard_pixelmapproducer + MapVecSize: 100000 + +} + +standard_cvnzlibmaker_east: +{ + module_type: ICARUSCVNZlibMaker + Module: "East" + ZLibFileName: "fwieler_pixel_maps" + InfoFileName: "fwieler_pixel_maps" + OutputDir: "./NewZlibMakerOutput" + Grid: "false" + PixelMapInput: "icaruspixelmapeast" + SetLog: false + ReverseViews: [false,false,false] + PlaneLimit: 500 + TDCLimit: 500 + verbose: true + UseSlice: true + SliceLabel: "pandoraGausCryoE" + TopologyHitsCut: 0 + GenieGenModuleLabel: "generator" + Verbose: true + UseBackTrackInfo: true + UseNuContainment: false + CosEfrac: 0.75 + NuEfrac: 0.25 + VolCut: 0. + SaveTree: false + PFParticleModuleLabel: "pandoraGausCryoE" + HitModuleLabel: "cluster3DCryoE" + T0Label: "pandoraGausCryoE" +} + +standard_cvnzlibmaker_west: +{ + module_type: ICARUSCVNZlibMaker + Module: West + ZLibFileName: "fwieler_pixel_maps" + InfoFileName: "fwieler_pixel_maps" + OutputDir: "./NewZlibMakerOutput" + PixelMapInput: "icaruspixelmapwest" + SetLog: false + ReverseViews: [false,false,false] + PlaneLimit: 500 + TDCLimit: 500 + verbose: true + UseSlice: true + SliceLabel: "pandoraGausCryoW" + TopologyHitsCut: 0 + GenieGenModuleLabel: "generator" + Verbose: true + UseBackTrackInfo: true + UseNuContainment: false + CosEfrac: 0.75 + NuEfrac: 0.25 + VolCut: 0. + SaveTree: false + PFParticleModuleLabel: "pandoraGausCryoW" + HitModuleLabel: "cluster3DCryoW" + T0Label: "pandoraGausCryoW" +} + + +standard_icarustfnethandler: +{ + LibPath: "" + #For the grid + #TFProtoBuf: "${INPUT_TAR_DIR_LOCAL}/icarus_cvn/srcs/icaruscode/icaruscode/ICARUSCVN/fcls/model/new_training_output1_flavour_40pixels_101layers_new_new_data_balanced_2-13-0.48.46.pb" + #Running Local + TFProtoBuf: "CVN/icarus_cvn_flavor.pb" + ChargeLogScale: false + NImageWires: 500 + NImageTDCs : 500 + ReverseViews: [false,false,false] + UseBundle: false + Inputs: ["view0", "view1", "view2"] + Outputs: ["flavour_1/Softmax"] + NInputs: 3 + NOutputs: 1 + verbose: true +} + +#IF RUNNING ON THE PIXEL MAPS NOT ON THE ROOT FILE +standard_icaruscvnevaluatorpm: +{ + module_type: ICARUSCVNEvaluator + Module: "West" + PixelMapModuleLabel: "" + PixelMapInput: "/pnfs/icarus/scratch/users/fwieler/pixelMaps_score2/W_set/" + SliceLabel: "" + PFParticleModuleLabel: "pandoraGausCryoW" + T0Label: "pandoraGausCryoW" + CVNType: "Tensorflow" + verbose: true + ICARUSTFHandler: { + tool_type: ICARUSTFNetHandler + @table::standard_icarustfnethandler + } + PixelMapProducer: @local::standard_pixelmapproducer +} + +#RUNNING ON THE ROOT FILE +standard_icaruscvnevaluatorslc_west: +{ + module_type: ICARUSCVNEvaluator + Module: "West" + PixelMapModuleLabel: "" + PixelMapInput: "" + SliceLabel: "pandoraGausCryoW" + PFParticleModuleLabel: "pandoraGausCryoW" + T0Label: "pandoraGausCryoW" + CVNType: "Tensorflow" + verbose: true + ICARUSTFHandler: { + tool_type: ICARUSTFNetHandler + @table::standard_icarustfnethandler + } + PixelMapProducer: @local::standard_pixelmapproducer +} + +standard_icaruscvnevaluatorslc_east: +{ + module_type: ICARUSCVNEvaluator + Module: "East" + PixelMapModuleLabel: "" + PixelMapInput: "" + SliceLabel: "pandoraGausCryoE" + PFParticleModuleLabel: "pandoraGausCryoE" + T0Label: "pandoraGausCryoE" + CVNType: "Tensorflow" + verbose: true + ICARUSTFHandler: { + tool_type: ICARUSTFNetHandler + @table::standard_icarustfnethandler + } + PixelMapProducer: @local::standard_pixelmapproducer +} + + +END_PROLOG diff --git a/icaruscode/ICARUSCVN/fcls/getTrueInfo.py b/icaruscode/ICARUSCVN/fcls/getTrueInfo.py new file mode 100644 index 000000000..31aedef73 --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/getTrueInfo.py @@ -0,0 +1,74 @@ +import re +import glob +import os + +def normalize(value): + if value > 2: + return 3 + return value + +def normalize2(value): + if value < 0: + return 1 + return 0 + +directory = '/pnfs/icarus/scratch/users/fwieler/pixelMaps_score2/W_set' # current directory +pattern = re.compile(r'e(\d+)') +files = list(glob.iglob(directory + "/**/*.gz", recursive=True)) + +output_lines = [] +extracted_e_values = [] +i = 0 + +for filename in files: + match = pattern.search(filename) + + if not match: + continue + i += 1 + print(f'Processing [{i}/{len(files)}]') + e_value = f"e{match.group(1)}" + extracted_e_values.append(e_value) + + file_ID = filename.split("/")[-1][:-3] + infofile = file_ID + '.info' + dir_path = os.path.dirname(filename) + + try: + with open(dir_path + "/" + infofile, 'r') as f: + info = f.readlines() + except FileNotFoundError: + continue + + if not len(info): + continue + + fInt = int(info[0].strip()) + flavour = fInt // 4 + if flavour == 3: + flavour = 2; + interaction = fInt % 4 + + fNuEnergy = float(info[1].strip()) + fLepEnergy = float(info[2].strip()) + fRecoNueEnergy = float(info[3].strip()) + fRecoNumuEnergy = float(info[4].strip()) + fEventWeight = float(info[6].strip()) + + fNuPDG = normalize2(int(info[7].strip())) + fNProton = normalize(int(info[8].strip())) + fNPion = normalize(int(info[9].strip())) + fNPizero = normalize(int(info[10].strip())) + fNNeutron = normalize(int(info[11].strip())) + + # Create a CSV line: eXX, flavour, interaction, energies, weights, normalized particles + line = f"{file_ID},{flavour},{interaction},{fNuEnergy},{fLepEnergy},{fRecoNueEnergy},{fRecoNumuEnergy},{fEventWeight},{fNuPDG},{fNProton},{fNPion},{fNPizero},{fNNeutron}" + output_lines.append(line) + +# Write all results to a file +with open("event_info_wset_score2.csv", "w") as outfile: + outfile.write("fid,flavour,interaction,fNuEnergy,fLepEnergy,fRecoNueEnergy,fRecoNumuEnergy,fEventWeight,fNuPDG,fNProton,fNPion,fNPizero,fNNeutron\n") + for line in output_lines: + outfile.write(line + "\n") + +print("\nSaved event info to 'event_info_wset_score2.csv'.") diff --git a/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNEvaluator.fcl b/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNEvaluator.fcl new file mode 100644 index 000000000..bcd5eb59f --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNEvaluator.fcl @@ -0,0 +1,62 @@ +#include "ICARUSCVNMapper.fcl" +#include "simulationservices_icarus.fcl" +#include "services_common_icarus.fcl" +#include "rootoutput_icarus.fcl" +#include "sam_icarus.fcl" +#include "services_icarus_simulation.fcl" +#include "signalservices_icarus.fcl" +#include "rootoutput_icarus.fcl" + +process_name: ICARUSCVNEVALUATOR + +services: +{ + TFileService: {} + RandomNumberGenerator: {} # required by fuzzyCluster + message: @local::icarus_message_services_prod #from services_common_icarus.fcl + FileCatalogMetadata: @local::art_file_catalog_mc #from sam_icarus.fcl + @table::icarus_common_services # include all ICARUS basic services + SignalShapingICARUSService: @local::icarus_signalshapingservice + ParticleInventoryService: @local::standard_particleinventoryservice + +} + +source: +{ + module_type: RootInput + maxEvents: -1 +} + + +physics: +{ + producers:{ + cvnCryoE: @local::standard_icaruscvnevaluatorslc_east + cvnCryoW: @local::standard_icaruscvnevaluatorslc_west + } + + filters: {} + + analyzers:{} + + reco: [ cvnCryoE, cvnCryoW ] + + stream: [ out ] + + trigger_paths: [ reco ] + + end_paths: [ stream ] +} + +outputs: +{ + out: + { + module_type: RootOutput + dataTier: "reconstructed" + SelectEvents: [ reco ] + fastCloning: true + fileName: @local::icarus_rootoutput.fileName + } +} + diff --git a/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNMapper.fcl b/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNMapper.fcl new file mode 100644 index 000000000..dec07750c --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/run_ICARUSCVNMapper.fcl @@ -0,0 +1,72 @@ +#include "ICARUSCVNMapper.fcl" +#include "simulationservices_icarus.fcl" +#include "services_common_icarus.fcl" +#include "rootoutput_icarus.fcl" +#include "services_icarus_simulation.fcl" +#include "sam_icarus.fcl" +#include "signalservices_icarus.fcl" + +process_name: ICARUSCVNMAPPER + +services: +{ + RandomNumberGenerator: {} # required by fuzzyCluster + message: @local::icarus_message_services_prod + FileCatalogMetadata: @local::art_file_catalog_mc # from sam_icarus.fcl + @table::icarus_basic_services + @table::icarus_common_services + SignalShapingICARUSService: @local::icarus_signalshapingservice + #SignalShapingICARUSService: @local::icarus_wirecalibration_minimum_services # from services_icarus_simulation.fcl + ParticleInventoryService: @local::standard_particleinventoryservice +} + +source: +{ + module_type: RootInput + maxEvents: -1 +} + +physics: +{ + producers:{ + icaruspixelmapeast: @local::standard_cvnmapper_east + icaruspixelmapwest: @local::standard_cvnmapper_west + } + + filters: {} + + analyzers:{} + + reco: [ icaruspixelmapeast, icaruspixelmapwest ] + + stream: [ out ] + + trigger_paths: [ reco ] + + end_paths: [ stream ] +} + +outputs: +{ + out: + { + module_type: RootOutput + dataTier: "reconstructed" + SelectEvents: [ reco ] + fileName: @local::icarus_rootoutput.fileName + } +} + +outputs.out.outputCommands: [ + "keep *_*_*_*", + "drop recob::Showers_*_*_*", + "drop recob::Tracks_*_*_*", + "drop recob::OpHits_*_*_*", + "drop recob::Clusters_*_*_*", + "drop recob::SpacePoints_*_*_*", + "drop anab::Calorimetrys_*_*_*", + "drop anab::anab::ParticleIDs_*_*_*", + "drop recob::OpFlashs_*_*_*" +] + + diff --git a/icaruscode/ICARUSCVN/fcls/run_ICARUSZlibMaker.fcl b/icaruscode/ICARUSCVN/fcls/run_ICARUSZlibMaker.fcl new file mode 100644 index 000000000..6d6d72c98 --- /dev/null +++ b/icaruscode/ICARUSCVN/fcls/run_ICARUSZlibMaker.fcl @@ -0,0 +1,45 @@ +#include "ICARUSCVNMapper.fcl" +#include "simulationservices_icarus.fcl" +#include "services_common_icarus.fcl" + +#include "sam_icarus.fcl" + +#include "services_icarus_simulation.fcl" +#include "signalservices_icarus.fcl" +#include "rootoutput_icarus.fcl" + +process_name: ZlibMaker + +services: +{ + TFileService: { fileName: "ICARUSZlibMaker.root" } + RandomNumberGenerator: {} # required by fuzzyCluster + message: @local::icarus_message_services_prod #from services_common_icarus.fcl + FileCatalogMetadata: @local::art_file_catalog_mc #from sam_icarus.fcl + @table::icarus_common_services # include all ICARUS basic services + SignalShapingICARUSService: @local::icarus_signalshapingservice + ParticleInventoryService: @local::standard_particleinventoryservice + BackTrackerService: @local::standard_backtrackerservice +} + +source: +{ + module_type: RootInput + maxEvents: -1 +} + +physics: +{ + producers:{} + + filters: {} + + analyzers:{ + mapeast:@local::standard_cvnzlibmaker_east + mapwest:@local::standard_cvnzlibmaker_west + } + + ana: [ mapeast, mapwest ] + end_paths: [ ana ] +} +services.BackTrackerService.BackTracker.SimChannelModuleLabel: "simdrift" diff --git a/icaruscode/ICARUSCVN/module_helpers/CMakeLists.txt b/icaruscode/ICARUSCVN/module_helpers/CMakeLists.txt new file mode 100644 index 000000000..e5df8cd56 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/CMakeLists.txt @@ -0,0 +1,43 @@ + +cet_make_library( + SOURCE + ICARUSPixelMapProducer.cxx + ICARUSICVNMapper.cxx + ICARUSICVNZlibMaker.cxx + ICARUSITFNetHandler.h + LIBRARY_NAME + icaruscode_icaruscvn_module_helpers + LIBRARIES + PUBLIC + larrecodnn::CVN_interfaces + larrecodnn::CVN_func + lardataobj::RecoBase + lardataobj::Simulation + art::Framework_Services_Registry + art::Framework_Core + art::Framework_Principal + art_root_io::tfile_support + art_root_io::TFileService_service + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + art_plugin_types::module + Boost::filesystem + art_root_io::tfile_support + nusimdata::SimulationBase + ROOT::Hist + lardata::DetectorPropertiesService + larcore::Geometry_Geometry_service + larcorealg::Geometry + + icaruscode::ICARUSCVN_tf + ${ZLIB_LIBRARIES} + PRIVATE + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service +) + +art_dictionary(DICTIONARY_LIBRARIES icaruscode_icaruscvn_module_helpers) + +install_headers() +install_source() diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.cxx b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.cxx new file mode 100644 index 000000000..e0b2477cb --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.cxx @@ -0,0 +1,142 @@ +#include "ICARUSICVNMapper.h" +#include "lardataobj/RecoBase/Slice.h" + +#include "canvas/Persistency/Common/FindManyP.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "lardata/Utilities/AssociationUtil.h" + +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/SubRun.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/View.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "art_root_io/TFileService.h" +#include "art_root_io/TFileDirectory.h" +#include "canvas/Persistency/Common/FindMany.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "fhiclcpp/ParameterSet.h" + +#include // std::min_element +#include // std::begin, std::end + +namespace lcvn +{ + template ICARUSICVNMapper::ICARUSICVNMapper(fhicl::ParameterSet const& pset) + : EDProducer{pset} + , fHitsModuleLabel(pset.get("HitsModuleLabel")) + , fClusterPMLabel(pset.get("ClusterPMLabel")) + , fMinClusterHits(pset.get("MinClusterHits")) + , fProducer(pset.get("PixelMapProducer")) + , fverbose(pset.get("verbose")) + , fUseSlice(pset.get("UseSlice")) + , fSliceLabel(pset.get("SliceLabel")) + , fPFParticleModuleLabel(pset.get("PFParticleModuleLabel")) + , fT0Label(pset.get("T0Label")) + , fMapVecSize(pset.get("MapVecSize")) + { + + produces>(fClusterPMLabel); + produces>(fClusterPMLabel); + } + + //-------------------------------------------------------------------------------------------------------------------------- + +template void ICARUSICVNMapper::produce(art::Event& evt) +{ + if(fverbose) std::cout << "============ Calling the function ICARUSICVNMapper::produce() ==============\n"; + + if(fUseSlice){ + if(fverbose) std::cout << "============ Calling the function ICARUSICVNMapper::produce() is using slices ==============\n"; + art::Handle< std::vector > SliceListHandle; + std::vector< art::Ptr > SliceList; + if( evt.getByLabel(fSliceLabel,SliceListHandle)) + art::fill_ptr_vector(SliceList,SliceListHandle); + + art::Handle< std::vector > PFPListHandle; + std::vector< art::Ptr > PFPList; + if( evt.getByLabel(fPFParticleModuleLabel,PFPListHandle)) + art::fill_ptr_vector(PFPList,PFPListHandle); + + art::FindManyP findManyHits(SliceListHandle, evt, fSliceLabel); + + art::FindManyP findManyPFPs(SliceListHandle, evt, fPFParticleModuleLabel); + art::FindManyP findManyT0s(PFPListHandle, evt, fT0Label); + + std::unique_ptr< std::vector > pmCol(new std::vector); + auto assn = std::make_unique< art::Assns >(); + + auto const detProp = art::ServiceHandle()->DataFor(evt); + + for(auto const& slice : SliceList){ + + if (fverbose) std::cout << "********* " << evt.run() << " " << evt.subRun() << " " << evt.id().event() << " " << slice->ID() << " **************\n"; + std::vector pfp_T0_vec; + if(findManyPFPs.isValid()){ + std::vector> slicePFPs = findManyPFPs.at(slice.key()); + if(slicePFPs.size()){ + for(auto const &pfp : slicePFPs){ + if(findManyT0s.isValid()){ + std::vector> T0_vec = findManyT0s.at(pfp.key()); + if(T0_vec.size()){ + for(auto const& T0 : T0_vec){ + pfp_T0_vec.push_back(T0->Time()); + } + } + } + } + } + } + + float min_T0 = 0.; + if(pfp_T0_vec.size()){ + min_T0 = *min_element(pfp_T0_vec.begin(), pfp_T0_vec.end()); + } + + if(findManyHits.isValid()){ + std::vector> slicehits = findManyHits.at(slice.key()); + fProducer.Set_fT0_value(min_T0); + PixelMap pm = fProducer.ICARUSCreateMap(detProp, slicehits); + auto nhits = fProducer.TotHits(); + pm.SetTotHits(nhits); + //pm.fSliceID = slice->ID(); + + if(nhits > fMinClusterHits && pmCol->size()push_back(pm); + util::CreateAssn(*this, evt, *pmCol, slice, *assn, fClusterPMLabel); + } + + } + } + //std::cout<size()<> hitlist; + auto hitListHandle = evt.getHandle>(fHitsModuleLabel); + if (hitListHandle) art::fill_ptr_vector(hitlist, hitListHandle); + + //Declaring containers for things to be stored in event + std::unique_ptr> pmCol(new std::vector); + + auto const detProp = art::ServiceHandle()->DataFor(evt); + PixelMap pm = fProducer.ICARUSCreateMap(detProp, hitlist); + auto nhits = fProducer.TotHits(); + pm.SetTotHits(nhits); + + if (nhits > fMinClusterHits) pmCol->push_back(pm); + + evt.put(std::move(pmCol), fClusterPMLabel); + } + } +} diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.h b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.h new file mode 100644 index 000000000..11c16ccf9 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.h @@ -0,0 +1,61 @@ +#ifndef LCVN_ICARUSICVNMAPPER_H +#define LCVN_ICARUSICVNMAPPER_H + +#include +#include +#include +#include +#include + +// Framework includes +#include "art/Framework/Core/EDProducer.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art_root_io/TFileDirectory.h" + +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Persistency/Common/Assns.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +// LArSoft includes +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/Simulation/SimChannel.h" + +#include "larrecodnn/CVN/func/PixelMap.h" +#include "larrecodnn/CVN/interfaces/PixelMapProducer.h" + +namespace lcvn { + + template + class ICARUSICVNMapper : public art::EDProducer { + public: + explicit ICARUSICVNMapper(fhicl::ParameterSet const& pset); + + void produce(art::Event& evt) override; + + protected: + /// Module lablel for input clusters + std::string fHitsModuleLabel; + + /// Instance lablel for cluster pixelmaps + std::string fClusterPMLabel; + + /// Minimum number of hits for cluster to be converted to pixel map + unsigned short fMinClusterHits; + + /// PixelMapProducer does the work for us + T fProducer; + + bool fverbose; + bool fUseSlice; + std::string fSliceLabel; + std::string fPFParticleModuleLabel; + std::string fT0Label; + unsigned int fMapVecSize; + }; + +} +#endif diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.cxx b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.cxx new file mode 100644 index 000000000..f7edb3040 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.cxx @@ -0,0 +1,52 @@ +#include "ICARUSICVNZlibMaker.h" +#include "lardataobj/RecoBase/Slice.h" + +#include "canvas/Persistency/Common/FindManyP.h" + +namespace fs = boost::filesystem; + +namespace lcvn { + + ICARUSICVNZlibMaker::ICARUSICVNZlibMaker(fhicl::ParameterSet const& pset) : EDAnalyzer(pset) + { + this->reconfigure(pset); + } + + //...................................................................... + ICARUSICVNZlibMaker::~ICARUSICVNZlibMaker() {} + + //...................................................................... + void ICARUSICVNZlibMaker::reconfigure(const fhicl::ParameterSet& pset) + { + fOutputDir = pset.get("OutputDir", ""); + fGrid = pset.get("Grid", false); + fModule = pset.get("Module"); + fPixelMapInput = pset.get("PixelMapInput"); + fSetLog = pset.get("SetLog"); + fReverseViews = pset.get>("ReverseViews"); + + fPlaneLimit = pset.get("PlaneLimit"); + fTDCLimit = pset.get("TDCLimit"); + fverbose = pset.get("verbose"); + fUseSlice = pset.get("UseSlice"); + fSliceLabel = pset.get("SliceLabel"); + } + + void ICARUSICVNZlibMaker::beginJob() + { + fImage = CVNImageUtils(fPlaneLimit, fTDCLimit, 3); + fImage.SetLogScale(fSetLog); + fImage.SetViewReversal(fReverseViews); + + if (fOutputDir != "") + out_dir = fOutputDir; + else + out_dir = "."; + + //Throw an error if the specified output directory doesn't exist + if (!fs::exists(out_dir)) + throw art::Exception(art::errors::FileOpenError) + << "Output directory " << out_dir << " does not exist!" << std::endl; + } + +} diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.h b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.h new file mode 100644 index 000000000..11753d6d5 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.h @@ -0,0 +1,87 @@ +#ifndef LCVN_ICARUSICVNZLIBMAKER_H +#define LCVN_ICARUSICVNZLIBMAKER_H + +#include +#include +#include +#include +#include +#include + +// C/C++ includes +#include +#include + +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/SubRun.h" +#include "art_root_io/TFileDirectory.h" +#include "boost/filesystem.hpp" + +// Framework includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/SubRun.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "canvas/Utilities/Exception.h" + +// Data products +#include "nusimdata/SimulationBase/GTruth.h" +#include "nusimdata/SimulationBase/MCFlux.h" +#include "nusimdata/SimulationBase/MCNeutrino.h" +#include "nusimdata/SimulationBase/MCTruth.h" +// #include "dunereco/FDSensOpt/FDSensOptData/EnergyRecoOutput.h" + +// CVN includes +#include "larrecodnn/CVN/func/AssignLabels.h" +#include "larrecodnn/CVN/func/CVNImageUtils.h" +#include "larrecodnn/CVN/func/InteractionType.h" +#include "larrecodnn/CVN/func/LArTrainingData.h" +#include "larrecodnn/CVN/func/PixelMap.h" + +// Compression +#include "math.h" +#include "zlib.h" + +#include "TH1.h" + +namespace lcvn { + + class ICARUSICVNZlibMaker : public art::EDAnalyzer { + public: + explicit ICARUSICVNZlibMaker(fhicl::ParameterSet const& pset); + ~ICARUSICVNZlibMaker(); + + void beginJob() override; + void analyze(const art::Event& evt) override {} + void reconfigure(const fhicl::ParameterSet& pset); + + protected: + std::string fOutputDir; + bool fGrid; + std::string fModule; + std::string fPixelMapInput; + bool fSetLog; + std::vector fReverseViews; + unsigned int fPlaneLimit; + unsigned int fTDCLimit; + bool fverbose; + bool fUseSlice; + std::string fSliceLabel; + + std::string out_dir; + CVNImageUtils fImage; + + template + void write_files(LArTrainingData td, std::string evtid) = delete; + + private: + std::string fZlibOutputName; + std::string fInfoOutputName; + }; +} +#endif diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h b/icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h new file mode 100644 index 000000000..27b715013 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h @@ -0,0 +1,29 @@ +///////////////////////////////////////////////////////////////////////////////////////////////// +/// \file ICARUSTFNetHandler.h +/// \brief ICARUSTFNetHandler for CVN +/// \author Varuna Meddage by looking at the format of larrecodnn/CVN/interfaces/ITFNetHandler.h +//////////////////////////////////////////////////////////////////////////////////////////////// + +#ifndef LCVN_ICARUSTFNETHANDLER_H +#define LCVN_ICARUSTFNETHANDLER_H + +#include +#include + +#include "larrecodnn/CVN/func/InteractionType.h" +#include "larrecodnn/CVN/func/PixelMap.h" + +namespace lcvn { + + /// Wrapper for caffe::Net which handles construction and prediction + class ICARUSITFNetHandler { + public: + virtual ~ICARUSITFNetHandler() noexcept = default; + /// Return prediction arrays for PixelMap + virtual std::vector> Predict(const PixelMap& pm, const int event, const std::string cryo) const = 0; + virtual std::vector> PredictFromArray(const std::vector &pa, const std::string event) const = 0; + }; + +} + +#endif // LCVN_ICARUSTFNETHANDLER_H diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.cxx b/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.cxx new file mode 100644 index 000000000..e5269e211 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.cxx @@ -0,0 +1,272 @@ +#include "ICARUSPixelMapProducer.h" + +namespace lcvn +{ + + template void ICARUSPixelMapProducer::ConvertLocaltoGlobal(geo::WireID wireid, unsigned int &globalWire, unsigned int &globalPlane) const + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::ConvertLocaltoGlobal() ==============\n"; + PixelMapProducer::ConvertLocaltoGlobal(wireid, globalWire, globalPlane); + } + + //------------------------------------------------------------------------------------------------------------------------------------------------------------------ + + template Boundary ICARUSPixelMapProducer::DefineBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T*>& cluster) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::DefineBoundary() ==============\n"; + + std::vector tmin_0; + std::vector tmin_1; + std::vector tmin_2; + + std::vector wire_0, bwire_0; + std::vector wire_1, bwire_1; + std::vector wire_2, bwire_2; + + std::vector tsum = {0., 0., 0.}; + std::vector tsize = {0., 0., 0.}; + + //std::cout << "===================== ICARUSPixelMapProducer::DefineBoundary() Number of hits : " << cluster.size() << "\n"; + + for (size_t iHit = 0; iHit < cluster.size(); ++iHit) { + U wraphit(*(cluster[iHit]), this->Threshold()); + Waveform wf = wraphit.GetWaveform(); + geo::WireID wireid = wraphit.GetID(); + + unsigned int tempWire = wireid.Wire; + unsigned int tempPlane = wireid.Plane; + + if (!this->MultipleDrifts()) ConvertLocaltoGlobal(wireid, tempWire, tempPlane); + + for (auto& pulse : wf) { + double min_tick = (double)INT_MAX; + for (auto& i : pulse) { + double temptdc = i.first; + if (this->MultipleDrifts()) ConvertLocaltoGlobalTDC(wireid, i.first, tempWire, tempPlane, temptdc); + if (temptdc < min_tick) min_tick = temptdc; + tsum[tempPlane] += temptdc; + tsize[tempPlane] += 1.; + } + + if (!(pulse.empty())) { + if (tempPlane == 0) { + tmin_0.push_back(min_tick); + wire_0.push_back(tempWire); + } + if (tempPlane == 1) { + tmin_1.push_back(min_tick); + wire_1.push_back(tempWire); + } + if (tempPlane == 2) { + tmin_2.push_back(min_tick); + wire_2.push_back(tempWire); + } + } + } // end loop over pulses on single wire + } // end loop over struck wires + + //std::cout << "===================== ICARUSPixelMapProducer::DefineBoundary() Reached end the hit loop ==================== " << cluster.size() << "\n"; + + double tmean_0 = tsum[0] / tsize[0]; + double tmean_1 = tsum[1] / tsize[1]; + double tmean_2 = tsum[2] / tsize[2]; + + for (int i = 0; i < (int)wire_0.size(); i++) { + if (std::abs(tmin_0[i] - tmean_0) < (double)this->TRes()) bwire_0.push_back(wire_0[i]); + } + for (int i = 0; i < (int)wire_1.size(); i++) { + if (std::abs(tmin_1[i] - tmean_1) < (double)this->TRes()) bwire_1.push_back(wire_1[i]); + } + for (int i = 0; i < (int)wire_2.size(); i++) { + if (std::abs(tmin_2[i] - tmean_2) < (double)this->TRes()) bwire_2.push_back(wire_2[i]); + } + + if (fverbose) std::cout << "Boundary wire vector sizes: " << bwire_0.size() << ", " << bwire_1.size() << ", " << bwire_2.size() << std::endl; + + int minwire_0 = 0; + int minwire_1 = 0; + int minwire_2 = 0; + auto minwireelement_0 = std::min_element(bwire_0.begin(), bwire_0.end()); + auto minwireelement_1 = std::min_element(bwire_1.begin(), bwire_1.end()); + auto minwireelement_2 = std::min_element(bwire_2.begin(), bwire_2.end()); + + if(fChangeWireNo){ + if (bwire_0.size() > 0) { + minwire_0 = *minwireelement_0 - 1; + if (fverbose) std::cout << "minwire 0: " << (*minwireelement_0 - 1) << std::endl; + } + if (bwire_1.size() > 0) { + minwire_1 = *minwireelement_1 - 1; + if (fverbose) std::cout << "minwire 1: " << (*minwireelement_1 - 1) << std::endl; + } + if (bwire_2.size() > 0) { + minwire_2 = *minwireelement_2 - 1; + if (fverbose) std::cout << "minwire 2: " << (*minwireelement_2 - 1) << std::endl; + } + } + + else{ + if (bwire_0.size() > 0) { + minwire_0 = *minwireelement_0; + if (fverbose) std::cout << "minwire 0: " << *minwireelement_0 << std::endl; + } + if (bwire_1.size() > 0) { + minwire_1 = *minwireelement_1; + if (fverbose) std::cout << "minwire 1: " << *minwireelement_1 << std::endl; + } + if (bwire_2.size() > 0) { + minwire_2 = *minwireelement_2; + if (fverbose) std::cout << "minwire 2: " << *minwireelement_2<< std::endl; + } + } + + this->SetTotHits(bwire_0.size() + bwire_1.size() + bwire_2.size()); + + Boundary bound(this->NWire(), this->TRes(), minwire_0, minwire_1, minwire_2, tmean_0, tmean_1, tmean_2); + + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::DefineBoundary() ==============\n"; + return bound; + } + + //------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template void ICARUSPixelMapProducer::ConvertLocaltoGlobalTDC(geo::WireID wireid, double localTDC, unsigned int &globalWire, unsigned int &globalPlane, double &globalTDC) const + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::ConvertLocaltoGlobalTDC() ==============\n"; + + // Fill hit table + size_t plane = wireid.Plane; + double time = localTDC; + if(wireid.TPC== 2 || wireid.TPC==3 ) + { + time = 6442.15 - localTDC; // 6500 = (2*960/0.4) + 850 the drif time plus the 850 tics addition + if(wireid.Plane==1) //the planes are changed in planes 1 and 2. The +60 ones are in the -60 ones of the other. + { + plane=2; + } + else if (wireid.Plane==2) + { + plane=1; + } + } + + if(wireid.TPC==1 || wireid.TPC==3) + { + if(plane==1 || plane == 2) + { + wireid.Wire = wireid.Wire + 2535; //2535 is the last part of the wires in cryos 0 an 2 fbefore the cut in z=0 + } + } + + globalWire = wireid.Wire; + globalPlane = plane; + globalTDC = time; + + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::ConvertLocaltoGlobalTDC() ==============\n"; + } + + //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template PixelMap ICARUSPixelMapProducer::ICARUSCreateMapGivenBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster, const Boundary& bound) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSCreateMapGivenBoundary() ==============\n"; + PixelMap pm(this->NWire(), this->NTdc(), bound); + + for(size_t iHit = 0; iHit < cluster.size(); ++iHit) + { + + U wraphit(*(cluster[iHit]), this->Threshold()); + Waveform wf = wraphit.GetWaveform(); + geo::WireID wireid = wraphit.GetID(); + + unsigned int tempWire = wireid.Wire; + unsigned int tempPlane = wireid.Plane; + + if(!this->MultipleDrifts()) + ConvertLocaltoGlobal(wireid, tempWire, tempPlane); + + for(auto &pulse: wf){ + // Leigh: Simple modification to unwrap the collection view wire plane + for(auto &i: pulse){ + const double pe = i.second; + double temptdc = i.first; + if(this->MultipleDrifts()) + ConvertLocaltoGlobalTDC(wireid, i.first, tempWire, tempPlane, temptdc); + + const unsigned int wire = tempWire; + const unsigned int wirePlane = tempPlane; + const double tdc = temptdc; + + pm.Add(wire, tdc, wirePlane, pe); + } + } + + } + pm.SetTotHits(this->TotHits()); + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::ICARUSCreateMapGivenBoundary() ==============\n"; + return pm; + } + + //-------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template PixelMap ICARUSPixelMapProducer::ICARUSCreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector>& cluster) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSCreateMap() [Art::Ptr form] ==============\n"; + std::vector newCluster; + for(const art::Ptr hit : cluster){ + newCluster.push_back(hit.get()); + } + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::ICARUSCreateMap() [Art::Ptr form] ==============\n"; + return ICARUSCreateMap(detProp, newCluster); + } + + //--------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template PixelMap ICARUSPixelMapProducer::ICARUSCreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSCreateMap() [Normal::Ptr form] ==============\n"; + Boundary bound = DefineBoundary(detProp, cluster); + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::ICARUSCreateMap() [Normal::Ptr form] ==============\n"; + return ICARUSCreateMapGivenBoundary(detProp, cluster, bound); + } + + //---------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template void ICARUSPixelMapProducer::Set_fT0_value(float value) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::Set_fT0_value() ==============\n"; + fT0 = value; + if(fverbose) std::cout << "============ Reached the end of the function ICARUSPixelMapProducer::Set_fT0_value() ==============\n"; + } + + //----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template PixelMap ICARUSPixelMapProducer::CreateMapGivenBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster, const Boundary& bound) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::CreateMapGivenBoundary() ==============\n"; + return PixelMapProducer::CreateMapGivenBoundary(detProp, cluster, bound); + } + + //------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------ + + template PixelMap ICARUSPixelMapProducer::CreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector>& cluster) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::CreateMap() ==============\n"; + return PixelMapProducer::CreateMap(detProp, cluster); + } + + //------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template PixelMap ICARUSPixelMapProducer::CreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster) + { + if(fverbose) std::cout << "============ Calling the function ICARUSPixelMapProducer::CreateMap() ==============\n"; + return PixelMapProducer::CreateMap(detProp, cluster); + } + + //------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------- + + template class ICARUSPixelMapProducer; + template class ICARUSPixelMapProducer; + template class ICARUSPixelMapProducer; + +} diff --git a/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h b/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h new file mode 100644 index 000000000..b12ad777c --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h @@ -0,0 +1,49 @@ +#ifndef LCVN_ICARUSPIXELMAPPRODUCER_H +#define LCVN_ICARUSPIXELMAPPRODUCER_H + +#include +#include +#include +#include +#include +#include + +#include "canvas/Persistency/Common/Ptr.h" +#include "larrecodnn/CVN/interfaces/PixelMapProducer.h" +#include "larrecodnn/CVN/func/PixelMap.h" + +namespace lcvn +{ + template class ICARUSPixelMapProducer : public PixelMapProducer + { + public: + ICARUSPixelMapProducer(unsigned int nWire, unsigned int nTdc, double tRes, double threshold = 0.):PixelMapProducer::PixelMapProducer(nWire, nTdc, tRes, threshold){std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSPixelMapProducer() ==============\n";} + ICARUSPixelMapProducer():PixelMapProducer::PixelMapProducer(){std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSPixelMapProducer() ==============\n";} + ICARUSPixelMapProducer(const fhicl::ParameterSet& pset):PixelMapProducer::PixelMapProducer(pset),fverbose(pset.get("verbose")),fChangeWireNo(pset.get("ChangeWireNo")),fReadoutSize(pset.get("ReadoutSize")),fShiftT(pset.get("ShiftT")),fInductionWires(pset.get("InductionWires")),fFlipInductionView(pset.get("FlipInductionView")),fUseT(pset.get("UseT")) {std::cout << "============ Calling the function ICARUSPixelMapProducer::ICARUSPixelMapProducer() ==============\n";} + Boundary DefineBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster) override; + void ConvertLocaltoGlobal(geo::WireID wireid, unsigned int &globalWire, unsigned int &globalPlane) const override; + void ConvertLocaltoGlobalTDC(geo::WireID wireid, double localTDC, unsigned int &globalWire, unsigned int &globalPlane, double &globalTDC) const override; + PixelMap CreateMapGivenBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster,const Boundary& bound) override; + PixelMap CreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector>& cluster) override; + PixelMap CreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster) override; + PixelMap ICARUSCreateMapGivenBoundary(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster,const Boundary& bound); + PixelMap ICARUSCreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector>& cluster); + PixelMap ICARUSCreateMap(detinfo::DetectorPropertiesData const& detProp,const std::vector< const T* >& cluster); + void Set_fT0_value(float value); + protected: + bool fverbose; + bool fChangeWireNo; + double fReadoutSize; // in time ticks + float fShiftT; // size of the back/front porch + int fInductionWires; // number of wires in the first induction plane + bool fFlipInductionView; // should we flip the induction view + bool fUseT; // + float fT0; // T0 coming from the PFP particles (in time ticks) + }; + + typedef ICARUSPixelMapProducer ICARUSPixelMapHitProducer; + typedef ICARUSPixelMapProducer ICARUSPixelMapWireProducer; + typedef ICARUSPixelMapProducer ICARUSPixelMapSimProducer; +} + +#endif // LCVN_ICARUSPIXELMAPPRODUCER_H diff --git a/icaruscode/ICARUSCVN/module_helpers/classes.h b/icaruscode/ICARUSCVN/module_helpers/classes.h new file mode 100644 index 000000000..9725b1cb6 --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/classes.h @@ -0,0 +1,8 @@ +#include "canvas/Persistency/Common/Wrapper.h" +#include "canvas/Persistency/Common/Assns.h" +#include "lardataobj/RecoBase/Slice.h" +#include "larrecodnn/CVN/func/PixelMap.h" +#include "larrecodnn/CVN/func/Result.h" +#include +#include +#include diff --git a/icaruscode/ICARUSCVN/module_helpers/classes_def.xml b/icaruscode/ICARUSCVN/module_helpers/classes_def.xml new file mode 100644 index 000000000..91995334f --- /dev/null +++ b/icaruscode/ICARUSCVN/module_helpers/classes_def.xml @@ -0,0 +1,23 @@ + + + + + + + + + + + + + + + + + + + + + + + diff --git a/icaruscode/ICARUSCVN/modules/CMakeLists.txt b/icaruscode/ICARUSCVN/modules/CMakeLists.txt new file mode 100644 index 000000000..0ecc5a12d --- /dev/null +++ b/icaruscode/ICARUSCVN/modules/CMakeLists.txt @@ -0,0 +1,100 @@ +cet_build_plugin(ICARUSCVNMapper art::EDProducer + LIBRARIES PRIVATE + larcorealg::Geometry + larcore::Geometry_Geometry_service + larsim::Simulation lardataobj::Simulation + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + lardata::Utilities + larevt::Filters + lardataobj::RawData + lardataobj::RecoBase + larreco::Calorimetry + larreco::RecoAlg + lardata::RecoObjects + larpandora::LArPandoraInterface + sbnobj::Common_CRT + nusimdata::SimulationBase + 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 canvas::canvas + art::Persistency_Provenance canvas::canvas + art::Utilities canvas::canvas + messagefacility::MF_MessageLogger + + fhiclcpp::fhiclcpp + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + ${ROOT_BASIC_LIB_LIST} + icaruscode_icaruscvn_module_helpers +) + +cet_build_plugin(ICARUSCVNZlibMaker art::EDAnalyzer + LIBRARIES PRIVATE + larcorealg::Geometry + larcore::Geometry_Geometry_service + larsim::Simulation lardataobj::Simulation + larsim::MCCheater_BackTrackerService_service + larsim::MCCheater_ParticleInventoryService_service + lardata::Utilities + larevt::Filters + lardataobj::RawData + lardataobj::RecoBase + larreco::Calorimetry + larreco::RecoAlg + lardata::RecoObjects + larpandora::LArPandoraInterface + sbnobj::Common_CRT + nusimdata::SimulationBase + 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 canvas::canvas + art::Persistency_Provenance canvas::canvas + art::Utilities canvas::canvas + messagefacility::MF_MessageLogger + + fhiclcpp::fhiclcpp + ROOT::Geom + ROOT::XMLIO + ROOT::Gdml + ${ROOT_BASIC_LIB_LIST} + icaruscode_icaruscvn_module_helpers +) + +cet_build_plugin(ICARUSTFNetHandler art::tool + LIBRARIES PRIVATE + larrecodnn::CVN_func + larrecodnn::CVN_interfaces + larrecodnn::ImagePatternAlgs_Tensorflow_TF + icaruscode::ICARUSCVN_tf + messagefacility::MF_MessageLogger + fhiclcpp::types + fhiclcpp::fhiclcpp + ROOT::Hist + icaruscode_icaruscvn_module_helpers + icaruscode::ICARUSCVN_tf +) + +cet_build_plugin(ICARUSCVNEvaluator art::EDProducer + LIBRARIES PRIVATE + larrecodnn::CVN_func + larrecodnn::CVN_interfaces + art_plugin_support::toolMaker + art::Framework_Core + art::Framework_Principal + art_root_io::TFileService_service + canvas::canvas + icaruscode_icaruscvn_module_helpers +) + +install_headers() +install_fhicl() +install_source() + diff --git a/icaruscode/ICARUSCVN/modules/ICARUSCVNEvaluator_module.cc b/icaruscode/ICARUSCVN/modules/ICARUSCVNEvaluator_module.cc new file mode 100644 index 000000000..70a84283f --- /dev/null +++ b/icaruscode/ICARUSCVN/modules/ICARUSCVNEvaluator_module.cc @@ -0,0 +1,338 @@ +//////////////////////////////////////////////////////////////////////// +// \file LArCVNEvaluator_module.cc +// \brief Producer module creating CVN neural net results +// \author Alexander Radovic - a.radovic@gmail.com +// Saul Alonso Monsalve - saul.alonso.monsalve@cern.ch +//////////////////////////////////////////////////////////////////////// + +// C/C++ includes +#include +#include +#include +#include +#include +#include +#include + +// Framework includes +#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/Services/Registry/ServiceHandle.h" +#include "art/Utilities/make_tool.h" +#include "art_root_io/TFileDirectory.h" +#include "art_root_io/TFileService.h" +#include "canvas/Persistency/Common/Assns.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h" +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h" +#include "lardata/Utilities/AssociationUtil.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/RecoBase/Slice.h" +#include "larrecodnn/CVN/func/AssignLabels.h" +#include "larrecodnn/CVN/func/InteractionType.h" +#include "larrecodnn/CVN/func/PixelMap.h" +#include "larrecodnn/CVN/func/Result.h" + +namespace lcvn { + +class ICARUSCVNEvaluator : public art::EDProducer { +public: + explicit ICARUSCVNEvaluator(fhicl::ParameterSet const &pset); + + void produce(art::Event &evt); + std::vector resolve_input_files(const std::string &pixelMapInpu); + std::vector load_from_folder(const std::string &folder_path); + std::vector decompress_pixelmap_gz(const std::string &filename); + void load_and_predict(const int event, + std::unique_ptr &handler, + const std::string &pixelMapInput); + +private: + // /// Module label for input pixel maps + std::string fModule; + std::string fPixelMapInput; + std::string fPixelMapModuleLabel; + art::InputTag fSliceLabel; + art::InputTag fPFParticleModuleLabel; + art::InputTag fT0Label; + std::string fCVNType; + bool fVerbose; + std::unique_ptr fTFHandler; + ICARUSPixelMapHitProducer fPMProducer; + // std::vector fEvaluators; +}; + +//....................................................................... +// ICARUSCVNEvaluator::ICARUSCVNEvaluator(fhicl::ParameterSet const &pset) +// : EDProducer{pset} { +// std::vector evalPsets = +// pset.get>("Evaluators"); + +// for (auto const &epset : evalPsets) { +// EvalConfig cfg{epset.get("Module"), +// epset.get("PixelMapInput"), +// epset.get("PixelMapModuleLabel"), +// epset.get("SliceLabel", ""), +// epset.get("PFParticleModuleLabel", ""), +// epset.get("T0Label", ""), +// epset.get("CVNType"), +// epset.get("verbose", false), +// art::make_tool( +// epset.get("ICARUSTFHandler")), +// ICARUSPixelMapHitProducer( +// epset.get("PixelMapProducer"))}; +// fEvaluators.push_back(std::move(cfg)); +// } +// produces>(); +// if (fEvaluators[0].fSliceLabel != "") { +// produces>(); +// } else { +// produces>(); +// } +// } +ICARUSCVNEvaluator::ICARUSCVNEvaluator(fhicl::ParameterSet const &pset) + : EDProducer{pset}, + fModule(pset.get("Module")), + fPixelMapInput(pset.get("PixelMapInput")), + fPixelMapModuleLabel(pset.get("PixelMapModuleLabel")), + fSliceLabel(pset.get("SliceLabel", "")), + fPFParticleModuleLabel(pset.get("PFParticleModuleLabel", "")), + fT0Label(pset.get("T0Label", "")), + fCVNType(pset.get("CVNType")), + fVerbose(pset.get("verbose", false)), + fTFHandler{art::make_tool(pset.get("ICARUSTFHandler"))}, + fPMProducer(pset.get("PixelMapProducer")) + { + produces>(); + if (fSliceLabel != "") { + produces>(); + } else { + produces>(); + } +} +//...................................................................... +void ICARUSCVNEvaluator::produce(art::Event &evt) { + + /// Define containers for the things we're going to produce + std::unique_ptr> resultCol(new std::vector); + auto assn1 = std::make_unique>(); + auto assn2 = std::make_unique>(); + bool isSlice = false; + // for (auto &cfg : fEvaluators) { + if (fCVNType == "TF" || fCVNType == "Tensorflow" || fCVNType == "TensorFlow") { + if (fSliceLabel != "") { // by default use slice information + if (fVerbose) std::cout << "Using slice information" << std::endl; + isSlice = true; + std::vector> slcList; + auto slcHandle = evt.getHandle>(fSliceLabel); + if (!slcHandle.isValid()) { + throw cet::exception("ICARUSCVNEvaluator") + << "Unable to get slices using label " << fSliceLabel; + } else { + if (fVerbose) std::cout << "Filling slice" << std::endl; + art::fill_ptr_vector(slcList, slcHandle); + } + art::Handle> PFPListHandle; + std::vector> PFPList; + if (evt.getByLabel(fPFParticleModuleLabel, PFPListHandle)) { + if (fVerbose) std::cout << "Filling pfp" << std::endl; + art::fill_ptr_vector(PFPList, PFPListHandle); + } + + art::FindManyP findManyHits(slcHandle, evt, + fSliceLabel); + art::FindManyP findManyPFPs( + slcHandle, evt, fPFParticleModuleLabel); + art::FindManyP findManyT0s(PFPListHandle, evt, fT0Label); + auto const detProp = + art::ServiceHandle() + ->DataFor(evt); + for (auto const &slice : slcList) { + if (fVerbose) + std::cout << "********* " << fModule << " ********* " + << evt.run() << " " << evt.subRun() << " " + << evt.id().event() << " " << slice->ID() + << " **************\n"; + std::vector pfp_T0_vec; + if (findManyPFPs.isValid()) { + std::vector> slicePFPs = findManyPFPs.at(slice.key()); + if (slicePFPs.size()) { + for (auto const &pfp : slicePFPs) { + if (findManyT0s.isValid()) { + std::vector> T0_vec = + findManyT0s.at(pfp.key()); + if (T0_vec.size()) { + for (auto const &T0 : T0_vec) { + pfp_T0_vec.push_back(T0->Time()); + } + } + } + } + } + } + + float min_T0 = 0.; + if (pfp_T0_vec.size()) { + min_T0 = *min_element(pfp_T0_vec.begin(), pfp_T0_vec.end()); + if (fVerbose) std::cout << "min t0: " << min_T0 << std::endl; + } + + if (findManyHits.isValid()) { + std::vector> slicehits = findManyHits.at(slice.key()); + if (fVerbose) std::cout << "slice key: " << slice.key() << " slice id: " << slice.id() << std::endl; + fPMProducer.Set_fT0_value(min_T0); + PixelMap pm = fPMProducer.ICARUSCreateMap(detProp, slicehits); + auto nhits = fPMProducer.TotHits(); + pm.SetTotHits(nhits); + // pm.fSliceID = slice->ID(); + std::string moduleSlcId = std::string(fModule) + " - slice id: " + std::to_string(slice.id()) + " - slice key: " + std::to_string(slice.key() + 1); + std::vector> output = fTFHandler->Predict(pm, evt.id().event(), moduleSlcId); + resultCol->emplace_back(output); + util::CreateAssn(*this, evt, *resultCol, slice, *assn1); + } else { + if (fVerbose) std::cout << "findManyHits: " << findManyHits.isValid() << std::endl; + } + } + } else { // Try to read pixel maps saved in the file + /// Load in the pixel maps + load_and_predict(evt.id().event(), fTFHandler, fPixelMapInput); + } + } else { + mf::LogError("ICARUSCVNEvaluator::produce") + << "CVN Type not in the allowed list: Tensorflow" << std::endl; + mf::LogError("ICARUSCVNEvaluator::produce") + << "Exiting without processing events" << std::endl; + return; + } + // } + for (const auto &r : *resultCol) { + for (const auto &row : r.fOutput) { + for (const auto &val : row) { + std::cout << val << " "; + } + std::cout << std::endl; + } + } + + evt.put(std::move(resultCol)); + if (isSlice) { + evt.put(std::move(assn1)); + } else { + evt.put(std::move(assn2)); + } +} + +// --- Main prediction logic --- +void ICARUSCVNEvaluator::load_and_predict( + const int event, std::unique_ptr &handler, + const std::string &pixelMapInput) { + auto files = resolve_input_files(pixelMapInput); + + int fileCount = 0; + for (const auto &f : files) { + ++fileCount; + std::cout << "Processing [" << fileCount << "/" << files.size() << "]" + << "\n"; + + try { + size_t lastSlash = f.find_last_of('/'); + std::string fileName = + (lastSlash != std::string::npos) ? f.substr(lastSlash + 1) : f; + + if (fileName.size() > 3 && + fileName.substr(fileName.size() - 3) == ".gz") { + fileName = fileName.substr(0, fileName.size() - 3); // Remove ".gz" + } + std::vector pixel_array = decompress_pixelmap_gz(f); + std::vector> output = + handler->PredictFromArray(pixel_array, fileName); + } catch (const std::exception &e) { + std::cerr << "Failed on " << f << ": " << e.what() << "\n"; + } + } +} + +std::vector +ICARUSCVNEvaluator::decompress_pixelmap_gz(const std::string &filename) { + constexpr size_t kPixelArraySize = + 3 * 500 * 500; // Set this to your actual image shape + std::vector pixel_array(kPixelArraySize); + + std::ifstream file(filename, std::ios::binary | std::ios::ate); + if (!file) + throw std::runtime_error("Cannot open file: " + filename); + auto size = file.tellg(); + file.seekg(0, std::ios::beg); + + std::vector compressed(size); + file.read(compressed.data(), size); + file.close(); + + uLongf dest_len = kPixelArraySize; + int res = uncompress(pixel_array.data(), &dest_len, + reinterpret_cast(compressed.data()), size); + + if (res != Z_OK || dest_len != kPixelArraySize) { + throw std::runtime_error( + "Decompression failed or pixel array size mismatch in " + filename); + } + + return pixel_array; +} + +// --- Load all .gz files from a folder --- +std::vector +ICARUSCVNEvaluator::load_from_folder(const std::string &folder_path) { + std::vector files; + for (const auto &entry : + std::filesystem::recursive_directory_iterator(folder_path)) { + if (entry.is_regular_file() && entry.path().extension() == ".gz") { + files.push_back(entry.path().string()); + } + } + if (files.empty()) + throw std::runtime_error("No .gz files found in folder: " + folder_path); + return files; +} + +// --- Load file list or single file --- +std::vector +ICARUSCVNEvaluator::resolve_input_files(const std::string &pixelMapInput) { + std::filesystem::path input(pixelMapInput); + + if (!std::filesystem::exists(input)) { + throw std::runtime_error("Input path does not exist: " + pixelMapInput); + } + + if (std::filesystem::is_regular_file(input)) { + if (input.extension() == ".gz") { + std::cout << "Single File" << std::endl; + return {pixelMapInput}; // Single .gz file + } else { + // Text file containing list of .gz files + std::ifstream listfile(pixelMapInput); + if (!listfile) + throw std::runtime_error("Cannot open list file: " + pixelMapInput); + std::vector paths; + std::string line; + while (std::getline(listfile, line)) { + if (!line.empty()) + paths.push_back(line); + } + return paths; + } + } else if (std::filesystem::is_directory(input)) { + return load_from_folder(pixelMapInput); // Folder + } + + throw std::runtime_error("Unsupported input type: " + pixelMapInput); +} + +DEFINE_ART_MODULE(lcvn::ICARUSCVNEvaluator) +} // namespace lcvn +//////////////////////////////////////////////////////////////////////// diff --git a/icaruscode/ICARUSCVN/modules/ICARUSCVNMapper_module.cc b/icaruscode/ICARUSCVN/modules/ICARUSCVNMapper_module.cc new file mode 100644 index 000000000..a14cd4f67 --- /dev/null +++ b/icaruscode/ICARUSCVN/modules/ICARUSCVNMapper_module.cc @@ -0,0 +1,17 @@ +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.cxx" +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSICVNMapper.h" + +#include "lardataobj/RecoBase/Hit.h" +#include "lardataobj/RecoBase/Wire.h" +#include "lardataobj/Simulation/SimChannel.h" + +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSPixelMapProducer.h" +#include "larrecodnn/CVN/func/PixelMap.h" + +namespace lcvn { + + typedef ICARUSICVNMapper ICARUSCVNMapper; + template class ICARUSICVNMapper; + +DEFINE_ART_MODULE(lcvn::ICARUSCVNMapper) +} diff --git a/icaruscode/ICARUSCVN/modules/ICARUSCVNZlibMaker_module.cc b/icaruscode/ICARUSCVN/modules/ICARUSCVNZlibMaker_module.cc new file mode 100644 index 000000000..5bb7129bf --- /dev/null +++ b/icaruscode/ICARUSCVN/modules/ICARUSCVNZlibMaker_module.cc @@ -0,0 +1,999 @@ +//////////////////////////////////////////////////////////////////////// +// \file CVNZlibMaker_module.cc +// \brief Analyzer module for creating CVN gzip file objects +// \author Jeremy Hewes - jhewes15@fnal.gov +// Saul Alonso-Monsalve - saul.alonso.monsalve@cern.ch +// - wrote the zlib code used in this module +//////////////////////////////////////////////////////////////////////// + +// C/C++ includes +#include +#include +#include + +#include "boost/filesystem.hpp" +#include "art/Framework/Principal/Event.h" +#include "art/Framework/Principal/Handle.h" +#include "art/Framework/Principal/SubRun.h" +#include "art_root_io/TFileDirectory.h" +#include "art_root_io/TFileService.h" + + +// Framework includes +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Framework/Principal/Event.h" +#include "canvas/Persistency/Common/Ptr.h" +#include "art/Framework/Principal/SubRun.h" +#include "canvas/Utilities/Exception.h" +#include "canvas/Persistency/Common/PtrVector.h" +#include "canvas/Persistency/Common/FindManyP.h" +#include "canvas/Persistency/Common/FindOneP.h" +#include "art/Framework/Services/Registry/ServiceHandle.h" + +// Data products +#include "nusimdata/SimulationBase/GTruth.h" +#include "nusimdata/SimulationBase/MCTruth.h" +#include "nusimdata/SimulationBase/MCFlux.h" +#include "nusimdata/SimulationBase/MCNeutrino.h" +#include "lardataobj/RecoBase/Slice.h" +#include "lardataobj/RecoBase/PFParticle.h" +#include "lardataobj/AnalysisBase/T0.h" +#include "lardataobj/RecoBase/Hit.h" +#include "larsim/MCCheater/ParticleInventoryService.h" +#include "lardata/DetectorInfoServices/DetectorClocksService.h" +#include "lardataalg/DetectorInfo/DetectorClocksData.h" +#include "larsim/MCCheater/BackTrackerService.h" +#include "lardataobj/Simulation/SimChannel.h" +#include "larpandora/LArPandoraInterface/LArPandoraHelper.h" +#include "lardataobj/RecoBase/PFParticleMetadata.h" +// #include "../../RecoUtils/RecoUtils.h" +#include "icaruscode/RecoUtils/RecoUtils.h" + +// #include "dunereco/FDSensOpt/FDSensOptData/EnergyRecoOutput.h" + +// CVN includes +#include "larrecodnn/CVN/func/AssignLabels.h" +#include "larrecodnn/CVN/func/LArTrainingData.h" +#include "larrecodnn/CVN/func/InteractionType.h" +#include "larrecodnn/CVN/func/PixelMap.h" +#include "larrecodnn/CVN/func/CVNImageUtils.h" +// #include "../module_helpers/ICARUSICVNZlibMaker.h" +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSICVNZlibMaker.h" + +// Compression +#include "zlib.h" +#include "math.h" + +#include "TH1.h" +#include "TTree.h" +#include "TFile.h" +#include "larcoreobj/SummaryData/POTSummary.h" + +namespace lcvn { + + class ICARUSCVNZlibMaker : public ICARUSICVNZlibMaker { + public: + + explicit ICARUSCVNZlibMaker(fhicl::ParameterSet const& pset); + ~ICARUSCVNZlibMaker(); + + void beginJob(); + void endSubRun(art::SubRun const &sr); + void analyze(const art::Event& evt); + void reconfigure(const fhicl::ParameterSet& pset); + + private: + + unsigned int fTopologyHitsCut; + std::string fGenieGenModuleLabel; + bool fVerbose; + bool fUseBackTrackInfo; + bool fUseNuContainment; + double fCosEfrac; + double fNuEfrac; + double fVolCut; + bool fSaveTree; + std::string fPFParticleModuleLabel; + std::string fHitModuleLabel; + std::string fT0Label; + + // ParticleInventoryService + art::ServiceHandle pi_serv; + + // BackTrackerService + art::ServiceHandle bt_serv; + + void write_files(LArTrainingNuData td, std::string evtid); + //bool Is_in_TPC(double stX, double stY, double stZ); + bool Is_in_Fiducial_Vol(double stX, double stY, double stZ); + void HitTruth(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, Int_t& truthid); + bool HitTruthId(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, Int_t& mcid); + bool TrackIdToMCTruth(Int_t const trkID, art::Ptr& mctruth); + double HitEfrmTrkID(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, Int_t truthid); + double HitTotE(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit); + bool Is_Slice_Nu(art::FindManyP const& assn, art::Ptr const& slce); + art::Ptr Get_Nu_like_PFP(art::FindManyP const& assn, art::Ptr const& slce); + double Get_Slice_Score(art::FindManyP const& assn, art::Ptr const& pfp); + std::vector Get_Best_Slice_ID_PFP_pdg(art::FindManyP const& pfp_assn, art::FindManyP const& metadata_assn, std::vector> const& slice_vec); + int Get_Slice_PFP_ID(art::FindManyP const& assn, art::Ptr const& slce); + int Get_Hit_Count(unsigned int tpc_no, unsigned int pln_no, std::vector> const& hit_vec); + double Get_Tot_nu_E(detinfo::DetectorClocksData const& clockData,art::Ptr const& mcneutrino, std::vector> const& hit_vec); + + double HitTotE_frm_given_orgin(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, std::string origin_name); + double HitTotE_frm_mctruth(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, art::Ptr const& my_mctruth); + double TotE_from_mctruth(detinfo::DetectorClocksData const& clockData, std::vector> const& hit_vec, art::Ptr const& my_mctruth); + + void Clear(); + + TH1D* hPOT; + double fPOT; + int fRun; + int fSubRun; + + TTree *fEventTree; + int frun; + int fsubrun; + int fevent; + int fsliceID; + double ftotsliceE; + double ftotsliceNuE; + double ftotsliceCosE; + double ftotsliceOthE; + double fsliceNuEfrac; + double fsliceCosEfrac; + double fsliceOthEfrac; + bool fIsSliceNu; + double fSliceScore; + bool fIsbestSlice; + int fbestpfppdg; + int fbestsliceID; + int fpfppdg; + bool fNuIsContained; + int fNuID; + int fNhits_tpc_0_pl_0; + int fNhits_tpc_0_pl_1; + int fNhits_tpc_0_pl_2; + int fNhits_tpc_1_pl_0; + int fNhits_tpc_1_pl_1; + int fNhits_tpc_1_pl_2; + int fNhits_total; + double fNuSlice_purity; + double fNuSlice_completeness; + double fT0; + std::string fZlibOutputName; + std::string fInfoOutputName; + }; + + //...................................................................... + ICARUSCVNZlibMaker::ICARUSCVNZlibMaker(fhicl::ParameterSet const& pset) + : ICARUSICVNZlibMaker(pset) + { + reconfigure(pset); + } + + //...................................................................... + ICARUSCVNZlibMaker::~ICARUSCVNZlibMaker() + { } + + //...................................................................... + void ICARUSCVNZlibMaker::reconfigure(const fhicl::ParameterSet& pset) + { + std::cout << "Reconfigure ZlibMaker" << std::endl; + ICARUSICVNZlibMaker::reconfigure(pset); + fTopologyHitsCut = pset.get("TopologyHitsCut"); + fGenieGenModuleLabel = pset.get("GenieGenModuleLabel"); + fVerbose = pset.get("Verbose"); + fUseBackTrackInfo = pset.get("UseBackTrackInfo"); + fUseNuContainment = pset.get("UseNuContainment"); + fCosEfrac = pset.get("CosEfrac"); + fNuEfrac = pset.get("NuEfrac"); + fVolCut = pset.get("VolCut"); + fSaveTree = pset.get("SaveTree"); + fPFParticleModuleLabel = pset.get("PFParticleModuleLabel"); + fHitModuleLabel = pset.get("HitModuleLabel"); + fT0Label = pset.get("T0Label"); + fZlibOutputName = pset.get("ZLibFileName"); + fInfoOutputName = pset.get("InfoFileName"); + } + + //...................................................................... + void ICARUSCVNZlibMaker::endSubRun(const art::SubRun & sr){ + std::cout << "EndSubRun" << std::endl; + std::string fPOTModuleLabel = "generator"; + fRun = sr.run(); + fSubRun = sr.subRun(); + + art::Handle< sumdata::POTSummary > potListHandle; + if(sr.getByLabel(fPOTModuleLabel,potListHandle)) + fPOT = potListHandle->totpot; + else + fPOT = 0.; + if(hPOT) hPOT->Fill(0.5, fPOT); + } + + //...................................................................... + + void ICARUSCVNZlibMaker::beginJob() + { + std::cout << "Begin Job ZlibMaker" << std::endl; + ICARUSICVNZlibMaker::beginJob(); + + art::ServiceHandle tfs; + hPOT = tfs->make("TotalPOT", "Total POT;; POT", 1, 0, 1); + } + + //...................................................................... + + void ICARUSCVNZlibMaker::analyze(const art::Event& evt) + { + std::cout << "Begin Analyze ZlibMaker" << std::endl; + + std::vector> pixelmaps; + art::InputTag itag1(fPixelMapInput, "cvnmap"); + auto h_pixelmaps = evt.getHandle>(itag1); + if (h_pixelmaps) + art::fill_ptr_vector(pixelmaps, h_pixelmaps); + + if (pixelmaps.size() == 0) return; + + // Get associated slice for each pixel map + art::FindOneP findOneSlice(h_pixelmaps, evt, itag1); + + AssignLabels labels; + TDNuInfo info; + double event_weight =1; + InteractionType interaction = kOther; + bool Nu_is_contained = false; + + std::vector> mctruth_list; + auto h_mctruth = evt.getHandle>(fGenieGenModuleLabel); + if (h_mctruth) + art::fill_ptr_vector(mctruth_list, h_mctruth); + + for(auto const& mctruth : mctruth_list){ + if(mctruth->Origin() == simb::kBeamNeutrino){ + simb::MCNeutrino true_neutrino = mctruth->GetNeutrino(); + + if(fUseNuContainment){ + if(Is_in_Fiducial_Vol(true_neutrino.Nu().Vx(),true_neutrino.Nu().Vy(),true_neutrino.Nu().Vz())){ + interaction = labels.GetInteractionType(true_neutrino); + labels.GetTopology(mctruth, fTopologyHitsCut); + float nu_energy = true_neutrino.Nu().E(); + float lep_energy = true_neutrino.Lepton().E(); + fNuID = true_neutrino.Nu().TrackId(); + info.SetTruthInfo(nu_energy, lep_energy, 0., event_weight); + info.SetTopologyInformation(labels.GetPDG(), labels.GetNProtons(), labels.GetNPions(), labels.GetNPizeros(), labels.GetNNeutrons(), labels.GetTopologyType(), labels.GetTopologyTypeAlt()); + Nu_is_contained = true; + break; + } + } + + else{ + interaction = labels.GetInteractionType(true_neutrino); + labels.GetTopology(mctruth, fTopologyHitsCut); + float nu_energy = true_neutrino.Nu().E(); + float lep_energy = true_neutrino.Lepton().E(); + fNuID = true_neutrino.Nu().TrackId(); + info.SetTruthInfo(nu_energy, lep_energy, 0., event_weight); + info.SetTopologyInformation(labels.GetPDG(), labels.GetNProtons(), labels.GetNPions(), labels.GetNPizeros(), labels.GetNNeutrons(), labels.GetTopologyType(), labels.GetTopologyTypeAlt()); + if(Is_in_Fiducial_Vol(true_neutrino.Nu().Vx(),true_neutrino.Nu().Vy(),true_neutrino.Nu().Vz())) Nu_is_contained = true; + break; + } + + } // select neutrino interactions + + + else if(mctruth->Origin() == simb::kCosmicRay){ + interaction = kCosmic; + break; + } // select cosmic interactions + } + + + /////////////////////// use slice section //////////////////////////////////////// + + if(fUseSlice){ + std::cout << "***************** Used slice method ****************\n"; + + art::Handle< std::vector > SliceListHandle; + std::vector< art::Ptr > SliceList; + if( evt.getByLabel(fSliceLabel,SliceListHandle)) + art::fill_ptr_vector(SliceList,SliceListHandle); + + art::Handle< std::vector > PFPListHandle; + std::vector< art::Ptr > PFPList; + if( evt.getByLabel(fPFParticleModuleLabel,PFPListHandle)) + art::fill_ptr_vector(PFPList,PFPListHandle); + + art::Handle< std::vector > HitListHandle; + std::vector< art::Ptr > HitList; + if( evt.getByLabel(fHitModuleLabel,HitListHandle)) + art::fill_ptr_vector(HitList,HitListHandle); + + + art::FindManyP findManyHits(SliceListHandle, evt, fSliceLabel); + art::FindManyP findManyPFPs(SliceListHandle, evt, fPFParticleModuleLabel); + art::FindManyP fm_pfpmd(PFPListHandle, evt, fPFParticleModuleLabel); + art::FindManyP findManyT0s(PFPListHandle, evt, fT0Label); + + if(fUseBackTrackInfo){ + std::cout << "***************** Used backtracker information to get neutrino information ****************\n"; + auto const clockData = art::ServiceHandle()->DataFor(evt); + + for(unsigned int i=0; iID() << "\n"; + AssignLabels labels; + TDNuInfo info; + double event_weight =1; + InteractionType interaction = kOther; + + fsliceID = slice->ID(); + if(findManyHits.isValid()){ + std::cout << "Found many hists at: " << slice.key() << std::endl; + std::vector> slice_hits = findManyHits.at(slice.key()); + double tot_slice_eng = 0; + double tot_slice_nu_eng = 0; + double tot_slice_cos_eng = 0; + + fNhits_tpc_0_pl_0 = Get_Hit_Count(0,0,slice_hits); + fNhits_tpc_0_pl_1 = Get_Hit_Count(0,1,slice_hits); + fNhits_tpc_0_pl_2 = Get_Hit_Count(0,2,slice_hits); + fNhits_tpc_1_pl_0 = Get_Hit_Count(1,0,slice_hits); + fNhits_tpc_1_pl_1 = Get_Hit_Count(1,1,slice_hits); + fNhits_tpc_1_pl_2 = Get_Hit_Count(1,2,slice_hits); + fNhits_total = slice_hits.size(); + + std::vector mc_truth_eng(mctruth_list.size(),0.); + for(auto const hit : slice_hits){ + //Int_t trkId; // newly deleted + tot_slice_eng += HitTotE(clockData,hit); + tot_slice_nu_eng += HitTotE_frm_given_orgin(clockData,hit,"nu"); // newly added + tot_slice_cos_eng += HitTotE_frm_given_orgin(clockData,hit,"cos"); // newly added + + // newly added section + + if(mctruth_list.size()){ + for(auto const truth : mctruth_list){ + if(truth->Origin() == simb::kBeamNeutrino){ + mc_truth_eng[truth.key()] += HitTotE_frm_mctruth(clockData,hit,truth); + } + } + } + } // loop over hits in the selected slice + + ftotsliceE = tot_slice_eng; + + if(tot_slice_eng > 0){ + std::cout << "Total energy : " << tot_slice_eng << "\n"; + std::cout << "Total cosmic energy : " << tot_slice_cos_eng << "\n"; + std::cout << "Total nu energy : " << tot_slice_nu_eng << "\n"; + std::cout << "Cosmic fraction : " << double(tot_slice_cos_eng)/tot_slice_eng << "\n"; + std::cout << "Neutrino fraction : " << double(tot_slice_nu_eng)/tot_slice_eng << "\n"; + + ftotsliceNuE = tot_slice_nu_eng; + ftotsliceCosE = tot_slice_cos_eng; + ftotsliceOthE = ftotsliceE - ftotsliceNuE - ftotsliceCosE; + fsliceNuEfrac = double(ftotsliceNuE)/ftotsliceE; + fsliceCosEfrac = double(ftotsliceCosE)/ftotsliceE; + fsliceOthEfrac = double(ftotsliceOthE)/ftotsliceE; + + ///////////////////////////////////////////////////// Check whether this has pandora T0 /////////////////////////////// + + std::vector pfp_T0_vec; + if(findManyPFPs.isValid()){ + std::cout << "0. Many PFPs is valid" << std::endl; + std::vector> slicePFPs = findManyPFPs.at(slice.key()); + if(slicePFPs.size()){ + for(auto const &pfp : slicePFPs){ + if(findManyT0s.isValid()){ + std::vector> T0_vec = findManyT0s.at(pfp.key()); + if(T0_vec.size()){ + for(auto const& T0 : T0_vec){ + pfp_T0_vec.push_back(T0->Time()); + } + } + } + } + } + } + + if(pfp_T0_vec.size()){ + fT0 = *min_element(pfp_T0_vec.begin(), pfp_T0_vec.end()); + } + + ////////////////////////////////////////////////////// //////////////////////////////////////////////////////////////// + + if((double(tot_slice_cos_eng)/tot_slice_eng) >= fCosEfrac){ + interaction = kCosmic; + } // select cosmic + else if((double(tot_slice_nu_eng)/tot_slice_eng) >= fNuEfrac){ + int index = -1; + double min_E = 0.; + std::cout << "Size of the truth energy list : " << mc_truth_eng.size() << "\n"; + for(unsigned int k=0; k min_E){ + //std::cout << "MC truth energy : " << mc_truth_eng[k] << "\n"; + art::Ptr mctruth = mctruth_list[k]; + simb::MCNeutrino true_neutrino = mctruth->GetNeutrino(); + //std::cout << "Neutrino vertext : " << true_neutrino.Nu().Vx() << " " << true_neutrino.Nu().Vy() << " " << true_neutrino.Nu().Vz() << "\n"; + + if(fUseNuContainment){ + if(Is_in_Fiducial_Vol(true_neutrino.Nu().Vx(),true_neutrino.Nu().Vy(),true_neutrino.Nu().Vz())){ + //std::cout << "================ Checking containment cut =================\n"; + index = k; + min_E = mc_truth_eng[k]; + fNuIsContained = true; + } + } + + else{ + //std::cout << "================ Not Checking containment cut =================\n"; + index = k; + min_E = mc_truth_eng[k]; + if(Is_in_Fiducial_Vol(true_neutrino.Nu().Vx(),true_neutrino.Nu().Vy(),true_neutrino.Nu().Vz())) fNuIsContained = true; + else fNuIsContained = false; + } + + /*if(fUseNuFullContainment){ + if(Is_in_Fiducial_Vol(true_neutrino.Nu().Vx(),true_neutrino.Nu().Vy(),true_neutrino.Nu().Vz())){ + index = k; + min_E = mc_truth_eng[k]; + } + }*/ + + + } + } + + if(index != -1){ + simb::MCNeutrino true_neutrino = mctruth_list[index]->GetNeutrino(); + interaction = labels.GetInteractionType(true_neutrino); + labels.GetTopology(mctruth_list[index], fTopologyHitsCut); + float nu_energy = true_neutrino.Nu().E(); + float lep_energy = true_neutrino.Lepton().E(); + fNuID = true_neutrino.Nu().TrackId(); + info.SetTruthInfo(nu_energy, lep_energy, 0., event_weight); + info.SetTopologyInformation(labels.GetPDG(), labels.GetNProtons(), labels.GetNPions(), labels.GetNPizeros(), labels.GetNNeutrons(), labels.GetTopologyType(), labels.GetTopologyTypeAlt()); + //fNuSlice_purity = double(Get_Tot_nu_E(clockData,mctruth_list[index],slice_hits))/tot_slice_eng; // newly deleted + fNuSlice_purity = double(mc_truth_eng[index])/tot_slice_eng; // newly added + //if(Get_Tot_nu_E(clockData,mctruth_list[index],HitList)>0)fNuSlice_completeness = double(Get_Tot_nu_E(clockData,mctruth_list[index],slice_hits))/Get_Tot_nu_E(clockData,mctruth_list[index],HitList); // newly deleted + if(TotE_from_mctruth(clockData, HitList, mctruth_list[index])) fNuSlice_completeness = double(mc_truth_eng[index])/TotE_from_mctruth(clockData, HitList, mctruth_list[index]); // newly added + std::cout << "Nu purity : " << fNuSlice_purity << " Nu completeness : " << fNuSlice_completeness << "\n"; + } + } // select nutrino + } // total slice energy is greater that 0 + } // valid hit association + + if(findManyPFPs.isValid()){ + std::cout << "1. Many PFPs is valid" << std::endl; + fIsSliceNu=Is_Slice_Nu(findManyPFPs,slice); + if(fIsSliceNu){ + if(fm_pfpmd.isValid()){ + fSliceScore = Get_Slice_Score(fm_pfpmd,Get_Nu_like_PFP(findManyPFPs,slice)); + if(slice->ID() == Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[0]) fIsbestSlice = true; + fbestpfppdg = Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[1]; + fbestsliceID = Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[0]; + fpfppdg = Get_Slice_PFP_ID(findManyPFPs, slice); + } + } + } + + LArTrainingNuData train(interaction, *pixelmaps[i], info); + std::string evtid = "r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"+std::to_string(evt.event())+"_sl"+std::to_string(slice->ID())+"_h"+std::to_string(time(0)); + std::cout << "0. Writing files" << std::endl; + write_files(train, evtid); + break; + //} // found the matching slice + } // find associated slice + } // loop over pixel maps + } // use backtrack inoformation + + + else{ + std::cout << "***************** Used truth level information to get neutrino information ****************\n"; + for(unsigned int i=0; i> slice_hits = findManyHits.at(slice.key()); + fNhits_tpc_0_pl_0 = Get_Hit_Count(0,0,slice_hits); + fNhits_tpc_0_pl_1 = Get_Hit_Count(0,1,slice_hits); + fNhits_tpc_0_pl_2 = Get_Hit_Count(0,2,slice_hits); + fNhits_tpc_1_pl_0 = Get_Hit_Count(1,0,slice_hits); + fNhits_tpc_1_pl_1 = Get_Hit_Count(1,1,slice_hits); + fNhits_tpc_1_pl_2 = Get_Hit_Count(1,2,slice_hits); + fNhits_total = slice_hits.size(); + } + + if(findManyPFPs.isValid()){ + fIsSliceNu=Is_Slice_Nu(findManyPFPs,slice); + + ////////////////////////////////////////// pandora T0 ////////////////////////////////////////////////////// + + std::vector pfp_T0_vec; + if(findManyPFPs.isValid()){ + std::vector> slicePFPs = findManyPFPs.at(slice.key()); + if(slicePFPs.size()){ + for(auto const &pfp : slicePFPs){ + if(findManyT0s.isValid()){ + std::vector> T0_vec = findManyT0s.at(pfp.key()); + if(T0_vec.size()){ + for(auto const& T0 : T0_vec){ + pfp_T0_vec.push_back(T0->Time()); + } + } + } + } + } + } + + if(pfp_T0_vec.size()){ + fT0 = *min_element(pfp_T0_vec.begin(), pfp_T0_vec.end()); + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////// + + if(fIsSliceNu){ + if(fm_pfpmd.isValid()){ + fSliceScore = Get_Slice_Score(fm_pfpmd,Get_Nu_like_PFP(findManyPFPs,slice)); + if(slice->ID() == Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[0]) fIsbestSlice = true; + fbestpfppdg = Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[1]; + fbestsliceID = Get_Best_Slice_ID_PFP_pdg(findManyPFPs,fm_pfpmd,SliceList)[0]; + fpfppdg = Get_Slice_PFP_ID(findManyPFPs, slice); + } + } + } + break; + //} + std::string evtid = "r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"+std::to_string(evt.event())+"_sl"+std::to_string(slice->ID())+"_h"+std::to_string(time(0)); + std::cout << "1. Writing files" << std::endl; + write_files(train, evtid); + } + + } + } // don't use slice information to extrac truth info + } // use slcie to make images + + else{ + std::cout << "***************** Used entire event ****************\n"; + LArTrainingNuData train(interaction, *pixelmaps[0], info); + // std::string evtid = "r"+std::to_string(evt.run())+"_s"+std::to_string(evt.subRun())+"_e"+std::to_string(evt.event())+"_h"+std::to_string(time(0)); + std::string evtid = "_e"+std::to_string(evt.event()); + std::cout << "2. Writing files" << std::endl; + write_files(train, evtid); + } // use event to make images + } + + //...................................................................... + + void ICARUSCVNZlibMaker::write_files(LArTrainingNuData td, std::string evtid) + { + std::string fileName = ""; + // cropped from 2880 x 500 to 500 x 500 here + std::vector pixel_array(3 * fPlaneLimit * fTDCLimit); + //fImage.DisableRegionSelection(); + // std::cout << "Wires: " << td.fPMap.NWire() << " tdc: " << td.fPMap.NTdc() << " pixel array size: " << pixel_array.size() << std::endl; + // std::cout << "n pixel: " << td.fPMap.NPixel() << "bound: " << td.fPMap.Bound() << " tothits: " << td.fPMap.GetTotHits() << std::endl; + + //USE THIS IF BREAKING ON MINMAXTDCS + //fImage.DisableRegionSelection(); + fImage.SetViewReversal(false, false, false); + fImage.SetPixelMapSize(td.fPMap.NWire(), td.fPMap.NTdc()); + fImage.ConvertPixelMapToPixelArray(td.fPMap, pixel_array); + + + ulong src_len = 3 * fPlaneLimit * fTDCLimit; // pixelArray length + ulong dest_len = compressBound(src_len); // calculate size of the compressed data + char* ostream = (char *) malloc(dest_len); // allocate memory for the compressed data + + int res = compress((Bytef *) ostream, &dest_len, (Bytef *) &pixel_array[0], src_len); + + // Buffer error + + if (res == Z_BUF_ERROR) + std::cout << "Buffer too small!" << std::endl; + + // Memory error + else if (res == Z_MEM_ERROR) + std::cout << "Not enough memory for compression!" << std::endl; + + // Compression ok + else { + // Create output files + // std::string image_file_name = "event_" + evtid + ".gz"; + // std::string info_file_name = "event_" + evtid + ".info"; + + std::random_device rd; // Non-deterministic seed (if available) + std::mt19937 gen(rd()); // Mersenne Twister engine + std::uniform_int_distribution<> dist1(0, 999); + std::uniform_int_distribution<> dist2(0, 1999999); + // Generate numbers + int r1 = dist1(gen); + int r2 = dist2(gen); + std::string dir = fGrid ? "" : "/"; + std::string image_file_name = fOutputDir + dir + fZlibOutputName + fileName + "_" + evtid + "_" + fModule + "_" + std::to_string(r1) + "_" + std::to_string(r2) + ".gz"; + std::string info_file_name = fOutputDir + dir + fInfoOutputName + fileName + "_" + evtid + "_" + fModule + "_" + std::to_string(r1) + "_" + std::to_string(r2) + ".info"; + + std::ofstream image_file (image_file_name, std::ofstream::binary); + std::ofstream info_file (info_file_name); + + std::cout << fZlibOutputName << " - file name: " << image_file_name << std::endl; + + if(image_file.is_open() && info_file.is_open()) { + + // Write compressed data to file + + image_file.write(ostream, dest_len); + + image_file.close(); // close file + + // Write records to file + + // Category + + info_file << td.fInt << std::endl; + //info_file << td.fInfo << std::endl; + info_file << td.fInfo; + info_file << td.fPMap.GetTotHits() << "," << fNhits_tpc_0_pl_0 << "," << fNhits_tpc_0_pl_1 << "," << fNhits_tpc_0_pl_2 << "," << fNhits_tpc_1_pl_0 << "," << fNhits_tpc_1_pl_1 << "," << fNhits_tpc_1_pl_2 << "," << fNhits_total << std::endl; + info_file << frun << "," << fsubrun << "," << fevent << "," << fsliceID << std::endl; + info_file << ftotsliceE << "," << ftotsliceNuE << "," << ftotsliceCosE << "," << ftotsliceOthE << "," << fsliceNuEfrac << "," << fsliceCosEfrac << "," << fsliceOthEfrac << "," << fNuSlice_purity << "," << fNuSlice_completeness << std::endl; + info_file << fIsSliceNu << "," << fSliceScore << "," << fIsbestSlice << "," << fbestpfppdg << "," << fbestsliceID << "," << fpfppdg << std::endl; + info_file << fNuIsContained << "," << fNuID << " " << fT0 << std::endl; + info_file.close(); // close file + } + else { + + if (image_file.is_open()) + image_file.close(); + else + throw art::Exception(art::errors::FileOpenError) + << "Unable to open file " << fZlibOutputName.c_str() << "!" << std::endl; + + if (info_file.is_open()) + info_file.close(); + else + throw art::Exception(art::errors::FileOpenError) + << "Unable to open file " << fInfoOutputName.c_str() << "!" << std::endl; + } + } + + free(ostream); // free allocated memory + + } // lcvn::CVNZlibMaker::write_files + + ////////////////////////////////////////////////////////////////////////////////////////////////////////// + + bool ICARUSCVNZlibMaker::Is_in_Fiducial_Vol(double stX, double stY, double stZ) + { + //if(fverbose) std::cout << "Is in Fiducial Vol Starting" << std::endl; + // double X_bound = 196.5 - fVolCut; + // double Y_bound = 200. - fVolCut; + // double Z_up_bound = 500. - fVolCut; + // double Z_low_bound = fVolCut; + + //USE ABS FOR X + double X_low_bound = 86.94; + double X_up_bound = 333.49; + double Y_low_bound = -156.86; + double Y_up_bound = 109.96; + double Z_low_bound = -864.95; + double Z_up_bound = 844.95; + + if((TMath::Abs(stX) >= X_low_bound && TMath::Abs(stX) <= X_up_bound) && + (stY >= Y_low_bound && stY <= Y_up_bound) && (stZ >=Z_low_bound && stZ <= Z_up_bound)) return true; + else return false; + } + + ///////////////////////////////////////////////////////////////////////////////////////////////////////// + + void ICARUSCVNZlibMaker::HitTruth( detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, Int_t& truthid){ + //if(fverbose) std::cout << "HitTruth Starting" << std::endl; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(clockData, hit); + if( !trackIDEs.size() ) return; + Float_t maxe = 0; + Int_t bestid = 0; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + if( trackIDEs[i].energy > maxe ) { + maxe = trackIDEs[i].energy; + bestid = trackIDEs[i].trackID; + } + } + truthid = bestid; + } + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////// + + bool ICARUSCVNZlibMaker::HitTruthId(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, Int_t& mcid){ + //if(fverbose) std::cout << "HitTruthId Starting" << std::endl; + mcid = std::numeric_limits::lowest(); + HitTruth(clockData,hit,mcid); + if( mcid > std::numeric_limits::lowest() ) return true; + else return false; + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////// + + bool ICARUSCVNZlibMaker::TrackIdToMCTruth(Int_t const trkID, art::Ptr& mctruth){ + bool matchFound = false; + //if(fverbose) std::cout << "TrackIdToMCTruth Starting" << std::endl; + try { + mctruth = pi_serv->TrackIdToMCTruth_P(trkID); + matchFound = true; + } catch(...) { + std::cout<<"Exception thrown matching TrackID "< const& hit, Int_t truthid){ + //if(fverbose) std::cout << "HitEfrmTrkID Starting" << std::endl; + double hit_E = 0.; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(clockData, hit); + if( !trackIDEs.size() ) return hit_E; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + if(trackIDEs[i].trackID == truthid) hit_E += trackIDEs[i].energy; + } + return hit_E; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::HitTotE(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit){ + //if(fverbose) std::cout << "HitTotE Starting" << std::endl; + double Tot_E = 0.; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(clockData, hit); + if( !trackIDEs.size() ) return Tot_E; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + Tot_E += trackIDEs[i].energy; + } + return Tot_E; + } + + ///////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + bool ICARUSCVNZlibMaker::Is_Slice_Nu(art::FindManyP const& assn, art::Ptr const& slce){ + //if(fverbose) std::cout << "Is_Slice_Nu Starting" << std::endl; + if(assn.isValid()){ + std::vector> slicePFPs = assn.at(slce.key()); + for(auto const &pfp : slicePFPs){ + if((pfp->PdgCode() == 12 || pfp->PdgCode() == 14)) return true; + } + } + return false; + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::Get_Slice_Score(art::FindManyP const& assn, art::Ptr const& pfp){ + //if(fverbose) std::cout << "Get_Slice_Score Starting" << std::endl; + double best_score = -9999; + if(assn.isValid() && pfp.isNonnull()){ + const std::vector> pfpMetaVec = assn.at(pfp->Self()); + if(pfpMetaVec.size()){ + for (auto const pfpMeta : pfpMetaVec){ + larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap(); + double score = propertiesMap.at("NuScore"); + if (score > best_score) best_score = score; + } + } + } + return best_score; + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + art::Ptr ICARUSCVNZlibMaker::Get_Nu_like_PFP(art::FindManyP const& assn, art::Ptr const& slce){ + art::Ptr result; + std::vector> slicePFPs = assn.at(slce.key()); + for(auto const &pfp : slicePFPs){ + if((pfp->PdgCode() == 12 || pfp->PdgCode() == 14)){ + result = pfp; + break; + } + } + return result; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + std::vector ICARUSCVNZlibMaker::Get_Best_Slice_ID_PFP_pdg(art::FindManyP const& pfp_assn, art::FindManyP const& metadata_assn, + std::vector> const& slice_vec){ + std::vector Sl_ID_pfp_pdf = {-999,-999}; + std::vector> nuPfpVec; + std::vector> nuSliceVec; + for(auto const &slice : slice_vec){ + std::vector> slicePFPs = pfp_assn.at(slice.key()); + for (auto const &pfp : slicePFPs){ + if ((pfp->PdgCode() == 12 || pfp->PdgCode() == 14)){ + nuPfpVec.push_back(pfp); + nuSliceVec.push_back(slice); + break; + } + } + } + + if(nuPfpVec.size()){ + float bestScore = -999.; + art::Ptr bestNuSlice; + art::Ptr bestNuPfp; + for(size_t i=0; i> pfpMetaVec = metadata_assn.at(nuPfpVec.at(i)->Self()); + for (auto const pfpMeta : pfpMetaVec) { + larpandoraobj::PFParticleMetadata::PropertiesMap propertiesMap = pfpMeta->GetPropertiesMap(); + score = propertiesMap.at("NuScore"); + if( score > bestScore ) { + bestScore = score; + bestNuPfp = nuPfpVec.at(i); + bestNuSlice = nuSliceVec.at(i); + } + } + } + + Sl_ID_pfp_pdf[0] = bestNuSlice->ID(); + Sl_ID_pfp_pdf[1] = bestNuPfp->PdgCode(); + } + return Sl_ID_pfp_pdf; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + int ICARUSCVNZlibMaker::Get_Slice_PFP_ID(art::FindManyP const& assn, art::Ptr const& slce){ + //if(fverbose) std::cout << "Get_Slice_PFP_ID Starting" << std::endl; + std::vector> slicePFPs = assn.at(slce.key()); + int result = -9999; + for(auto const &pfp : slicePFPs){ + if((pfp->PdgCode() == 12 || pfp->PdgCode() == 14)){ + result = pfp->PdgCode(); + break; + } + } + return result; + } + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + int ICARUSCVNZlibMaker::Get_Hit_Count(unsigned int tpc_no, unsigned int pln_no, std::vector> const& hit_vec){ + //if(fverbose) std::cout << "Get_Hit_Count Starting" << std::endl; + int n_hits = 0; + if(hit_vec.size()){ + for(auto const hit : hit_vec){ + if(hit->WireID().TPC == tpc_no){ + if(hit->WireID().Plane == pln_no){ + n_hits++; + } + } + } + } + + return n_hits; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::Get_Tot_nu_E(detinfo::DetectorClocksData const& clockData,art::Ptr const& mcneutrino, std::vector> const& hit_vec){ + double tot_nu_E = 0; + for(auto const hit : hit_vec){ + Int_t trkId; + if(HitTruthId(clockData,hit,trkId)){ + art::Ptr mctruth; + if(TrackIdToMCTruth(trkId,mctruth)){ + if(mctruth->Origin() == simb::kBeamNeutrino){ + if(mctruth == mcneutrino){ + tot_nu_E += HitEfrmTrkID(clockData,hit,trkId); + } + } + } + } + } + return tot_nu_E; + } + + + ////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::HitTotE_frm_given_orgin(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, std::string origin_name){ + double Tot_E = 0.; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(clockData, hit); + if( !trackIDEs.size() ) return Tot_E; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + art::Ptr mctruth; + if(TrackIdToMCTruth(trackIDEs[i].trackID,mctruth)){ + if(origin_name == "nu"){ + if(mctruth->Origin() == simb::kBeamNeutrino) Tot_E += trackIDEs[i].energy; + } + else if(origin_name == "cos"){ + if(mctruth->Origin() == simb::kCosmicRay) Tot_E += trackIDEs[i].energy; + } + } + + } + return Tot_E; + } + + //////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::HitTotE_frm_mctruth(detinfo::DetectorClocksData const& clockData, art::Ptr const& hit, art::Ptr const& my_mctruth){ + double Tot_E = 0.; + std::vector trackIDEs = bt_serv->HitToTrackIDEs(clockData, hit); + if( !trackIDEs.size() ) return Tot_E; + for(size_t i = 0; i < trackIDEs.size(); ++i){ + art::Ptr mctruth; + if(TrackIdToMCTruth(trackIDEs[i].trackID,mctruth)){ + if(my_mctruth.isNonnull()){ + if(mctruth == my_mctruth){ + Tot_E += trackIDEs[i].energy; + } + } + } + } + return Tot_E; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + double ICARUSCVNZlibMaker::TotE_from_mctruth(detinfo::DetectorClocksData const& clockData, std::vector> const& hit_vec, art::Ptr const& my_mctruth){ + double tot_E = 0; + if(hit_vec.size() && my_mctruth.isNonnull()){ + for(auto const hit : hit_vec){ + tot_E += HitTotE_frm_mctruth(clockData,hit,my_mctruth); + } + } + return tot_E; + } + + /////////////////////////////////////////////////////////////////////////////////////////////////////////////// + + void ICARUSCVNZlibMaker::Clear(){ + if(fverbose) std::cout << "Clear ZlibMaker " << std::endl; + frun = -9999; + fsubrun = -9999; + fevent = -9999; + fsliceID = -9999; + ftotsliceE = -9999; + ftotsliceNuE = -9999; + ftotsliceCosE = -9999; + ftotsliceOthE = -9999; + fsliceNuEfrac = -9999; + fsliceCosEfrac = -9999; + fsliceOthEfrac = -9999; + fIsSliceNu = false; + fSliceScore = -9999; + fIsbestSlice = false; + fbestpfppdg = -9999; + fbestsliceID = -9999; + fpfppdg = -9999; + fNuIsContained = false; + fNuID = -9999; + fNhits_tpc_0_pl_0 = -9999; + fNhits_tpc_0_pl_1 = -9999; + fNhits_tpc_0_pl_2 = -9999; + fNhits_tpc_1_pl_0 = -9999; + fNhits_tpc_1_pl_1 = -9999; + fNhits_tpc_1_pl_2 = -9999; + fNhits_total = -9999; + fNuSlice_purity = -9999.; + fNuSlice_completeness = -9999.; + fT0 = 0.; + } + + + DEFINE_ART_MODULE(ICARUSCVNZlibMaker) +} // namespace cvn diff --git a/icaruscode/ICARUSCVN/modules/ICARUSTFNetHandler_tool.cc b/icaruscode/ICARUSCVN/modules/ICARUSTFNetHandler_tool.cc new file mode 100644 index 000000000..0538a4e08 --- /dev/null +++ b/icaruscode/ICARUSCVN/modules/ICARUSTFNetHandler_tool.cc @@ -0,0 +1,298 @@ +//////////////////////////////////////////////////////////////////////////////////// +/// \file ICARUSTFNetHandler.cxx +/// \brief ICARUSTFNetHandler for CVN +/// \author Varuna Meddage (copied from +/// larrecodnn/CVN/tools/TFNetHandler_tool.cc) modified by Felipe Wieler for +/// Icarus CVN +/////////////////////////////////////////////////////////////////////////////////// + +#include "cetlib/getenv.h" +#include +#include +#include // for std::getenv +#include +#include + +#include "art/Utilities/ToolMacros.h" +#include "canvas/Utilities/Exception.h" +#include "fhiclcpp/ParameterSet.h" +#include "messagefacility/MessageLogger/MessageLogger.h" + +#include "icaruscode/ICARUSCVN/module_helpers/ICARUSITFNetHandler.h" +#include "icaruscode/ICARUSCVN/tf/tf_graph.h" +#include "larrecodnn/CVN/func/CVNImageUtils.h" + +namespace lcvn { + +/// Wrapper for caffe::Net which handles construction and prediction +class ICARUSTFNetHandler : public ICARUSITFNetHandler { +public: + /// Constructor which takes a pset with DeployProto and ModelFile fields + explicit ICARUSTFNetHandler(const fhicl::ParameterSet &pset); + + /// Return prediction arrays for PixelMap + std::vector> Predict(const PixelMap &pm, const int event, const std::string cryo) const override; + std::vector> PredictFromArray(const std::vector &pa, const std::string event) const override; + +private: + std::string fLibPath; ///< Library path (typically dune_pardata...) + std::string + fTFProtoBuf; ///< location of the tf .pb file in the above path or the + ///< directory containing model files in SavedModel format + ///< (set UseBundle = true in this case) + bool fUseLogChargeScale; ///< Is the charge using a log scale? + unsigned int fImageWires; ///< Number of wires for the network to classify + unsigned int fImageTDCs; ///< Number of tdcs for the network to classify + std::vector fReverseViews; ///< Do we need to reverse any views? + bool fUseBundle; ///< Use a bundled model saved in the SavedModel format from + ///< Tensorflow + std::vector fInputs; + std::vector fOutputs; + int fNInputs; + int fNOutputs; + std::unique_ptr fTFGraph; ///< Tensorflow graph + bool fverbose; +}; + +std::string expandEnvVars(const std::string& text) { + static std::regex env("\\$\\{([^}]+)\\}"); + std::smatch match; + std::string result = text; + while (std::regex_search(result, match, env)) { + const char* val = std::getenv(match[1].str().c_str()); + result.replace(match.position(0), match.length(0), + val ? val : ("${" + match[1].str() + "}")); + } + std::cout << result << std::endl; + return result; +} + +ICARUSTFNetHandler::ICARUSTFNetHandler(const fhicl::ParameterSet &pset) + // : fLibPath(cet::getenv(pset.get("LibPath", ""))) + : fTFProtoBuf(expandEnvVars(pset.get("TFProtoBuf"))), + fUseLogChargeScale(pset.get("ChargeLogScale")), + fImageWires(pset.get("NImageWires")), + fImageTDCs(pset.get("NImageTDCs")), + fReverseViews(pset.get>("ReverseViews")), + fUseBundle(pset.get("UseBundle")), + fInputs(pset.get>("Inputs")), + fOutputs(pset.get>("Outputs")), + fNInputs(pset.get("NInputs")), fNOutputs(pset.get("NOutputs")), + fverbose(pset.get("verbose")) { + + // Construct the TF Graph object. The empty vector {} is used since the + // protobuf file gives the names of the output layer nodes + cet::search_path sp("FW_SEARCH_PATH"); // read environment variable + std::string TFProtoBuf = sp.find_file(fTFProtoBuf.c_str()); // relative path + + mf::LogInfo("TFNetHandler") + << "Loading network: " << TFProtoBuf << std::endl; + fTFGraph = tf::Graph::create(TFProtoBuf.c_str(), fInputs, fOutputs, + fUseBundle, fNInputs, fNOutputs); + if (!fTFGraph) { + art::Exception(art::errors::Unknown) + << "Tensorflow model not found or incorrect"; + } +} + +// Check the network outputs +bool check(const std::vector> &outputs) { + if (outputs.size() == 1) + return true; + size_t aux = 0; + for (size_t o = 0; o < outputs.size(); ++o) { + size_t aux2 = 0; + + for (size_t i = 0; i < outputs[o].size(); ++i) + if (outputs[o][i] == 0.0 || outputs[o][i] == 1.0) + aux2++; + if (aux2 == outputs[o].size()) + aux++; + } + return aux == outputs.size() ? false : true; +} + +// Fill outputs with value -3 +void fillEmpty(std::vector> &outputs) { + for (auto &output : outputs) { + output.assign(output.size(), -3.0); + } + + return; +} + +std::vector> +ICARUSTFNetHandler::Predict(const PixelMap &pm, const int event, const std::string cryo) const { + + CVNImageUtils imageUtils(fImageWires, fImageTDCs, 3); + // Configure the image utility + imageUtils.SetViewReversal(fReverseViews); + imageUtils.SetImageSize(fImageWires, fImageTDCs, 3); + imageUtils.SetLogScale(fUseLogChargeScale); + imageUtils.SetPixelMapSize(pm.NWire(), pm.NTdc()); + + ImageVectorF thisImage; + imageUtils.ConvertPixelMapToImageVectorF(pm, thisImage); + std::vector vecForTF; + + vecForTF.push_back(thisImage); + std::vector> target_names; + if (fNOutputs == 1){ + target_names = { + {"CC Numu", "CC Nue", "Cosmic", "NC"}}; + } + else if (fNOutputs == 3){ + target_names = { + {"neutrino", "antineutrino"}, + {"CC Numu", "CC Nue", "NC"}, + {"CC QE", "CC Res", "CC DIS", "CC Other", "NC"}}; + } + else if(fNOutputs == 7){ + target_names = { + {"neutrino", "antineutrino", "NULL"}, + {"CC Numu", "CC Nue", "CC Nutau", "NC"}, + {"CC QE", "CC Res", "CC DIS", "CC Other", "NULL"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}}; + } + + std::vector>> + cvnResults; // shape(samples, #outputs, output_size) + bool status = false; + + int counter = 0; + + while (status == false) { // do until it gets a correct result + cvnResults = fTFGraph->run(vecForTF); + status = check(cvnResults[0]); + + counter++; + if (counter == 10) { + std::cout << "Error, CVN never outputing a correct result. Filling " + "result with zeros."; + std::cout << std::endl; + fillEmpty(cvnResults[0]); + break; + } + }; + + //TO MAINTAIN THE SBND FORMAT FOR THE CAFMAKER FILES. + for (auto &outer : cvnResults) { + for (auto &inner : outer) { + inner.insert(inner.begin() + 2, 0.0f); + } + } + + if (fverbose) { + + std::string s = "Classifier summary: e" + std::to_string(event); + int output_index = 0; + for (auto const &output : cvnResults[0]) { + s += "\nOutput " + std::to_string(output_index) + ": "; + for (auto const v : output) + s += std::to_string(v) + ", "; + s += cryo; + std::cout << s << std::endl; + output_index++; + } + + std::ofstream outfile("classifier_summary.txt", std::ios::app); // Save to file + + outfile << s; + outfile << std::endl; + + outfile.close(); + } + + return cvnResults[0]; +} + +std::vector> +ICARUSTFNetHandler::PredictFromArray(const std::vector &pa, const std::string event) const { + + CVNImageUtils imageUtils(fImageWires, fImageTDCs, 3); + // Configure the image utility + imageUtils.SetViewReversal(fReverseViews); + imageUtils.SetImageSize(fImageWires, fImageTDCs, 3); + imageUtils.SetLogScale(fUseLogChargeScale); + + ImageVectorF thisImage; + imageUtils.ConvertPixelArrayToImageVectorF(pa, thisImage); + std::vector vecForTF; + + vecForTF.push_back(thisImage); + std::vector> target_names; + if (fNOutputs == 1){ + target_names = { + {"CC Numu", "CC Nue", "NC"}}; + } + else if (fNOutputs == 3){ + target_names = { + {"neutrino", "antineutrino"}, + {"CC Numu", "CC Nue", "NC"}, + {"CC QE", "CC Res", "CC DIS", "CC Other", "NC"}}; + } + else if(fNOutputs == 7){ + target_names = { + {"neutrino", "antineutrino", "NULL"}, + {"CC Numu", "CC Nue", "CC Nutau", "NC"}, + {"CC QE", "CC Res", "CC DIS", "CC Other", "NULL"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}, + {"0", "1", "2", ">2"}}; + } + + std::vector>> + cvnResults; // shape(samples, #outputs, output_size) + bool status = false; + + int counter = 0; + + while (status == false) { // do until it gets a correct result + cvnResults = fTFGraph->run(vecForTF); + status = check(cvnResults[0]); + + counter++; + if (counter == 10) { + std::cout << "Error, CVN never outputing a correct result. Filling " + "result with zeros."; + std::cout << std::endl; + fillEmpty(cvnResults[0]); + break; + } + }; + + if (fverbose) { + + std::string s = "Classifier summary: " + event; + int output_index = 0; + for (auto const &output : cvnResults[0]) { + s += "\nOutput " + std::to_string(output_index) + ": "; + auto max_it = std::max_element(output.begin(), output.end()); + int argmax = std::distance(output.begin(), max_it); + if (output_index == 0) + argmax = std::round(output[0]); + for (auto const v : output) + s += std::to_string(v) + ", "; + s += " " + target_names[output_index][argmax] + " - " + std::to_string(*max_it); + std::cout << event << target_names[output_index][argmax] << " - " << *max_it + << std::endl; + output_index++; + } + + std::ofstream outfile("classifier_summary.txt", std::ios::app); // Save to file + + outfile << s; + outfile << std::endl; + + outfile.close(); + } + + return cvnResults[0]; +} + +} // namespace lcvn +DEFINE_ART_CLASS_TOOL(lcvn::ICARUSTFNetHandler) diff --git a/icaruscode/ICARUSCVN/tf/CMakeLists.txt b/icaruscode/ICARUSCVN/tf/CMakeLists.txt new file mode 100644 index 000000000..6408286fc --- /dev/null +++ b/icaruscode/ICARUSCVN/tf/CMakeLists.txt @@ -0,0 +1,9 @@ +cet_make_library(SOURCE + tf_graph.cc + LIBRARIES PRIVATE + TensorFlow::cc + TensorFlow::framework +) + +install_headers() +install_source() diff --git a/icaruscode/ICARUSCVN/tf/tf_graph.cc b/icaruscode/ICARUSCVN/tf/tf_graph.cc new file mode 100644 index 000000000..77d87997b --- /dev/null +++ b/icaruscode/ICARUSCVN/tf/tf_graph.cc @@ -0,0 +1,423 @@ +//////////////////////////////////////////////////////////////////////////////////////////////////// +// Class: Graph +// Authors: R.Sulej (Robert.Sulej@cern.ch), from DUNE, FNAL/NCBJ, Sept. 2017 +// P.Plonski, from DUNE, WUT, Sept. 2017 +// T.Cai (tejinc@yorku.ca) from DUNE, YorkU, March 2022 +// B.N.Nayak (nayakb@uci.edu) from DUNE, BNL, November 2022 +// +// Iterface to run Tensorflow graph saved to a file. First attempts, quite functional. +// +//////////////////////////////////////////////////////////////////////////////////////////////////// + +#include "tf_graph.h" + +#include "tensorflow/cc/saved_model/loader.h" +#include "tensorflow/cc/saved_model/tag_constants.h" +#include "tensorflow/core/platform/env.h" + +#include "tensorflow/core/public/session_options.h" + +// ------------------------------------------------------------------- +tf::Graph::Graph(const char* graph_file_name, + const std::vector& inputs, + const std::vector& outputs, + bool& success, + bool use_bundle, + int ninputs, + int noutputs) +{ + fUseBundle = use_bundle; + success = false; // until all is done correctly + + n_inputs = ninputs; + n_outputs = noutputs; + + // Force tf to only use a single core so it doesn't eat batch farms + tensorflow::SessionOptions options; + tensorflow::ConfigProto& config = options.config; + config.set_inter_op_parallelism_threads(1); + config.set_intra_op_parallelism_threads(1); + config.set_use_per_session_threads(false); + + auto status = tensorflow::NewSession(options, &fSession); + if (!status.ok()) { + std::cout << status.ToString() << std::endl; + return; + } + + if (fUseBundle) { + + fBundle = new tensorflow::SavedModelBundle(); + status = tensorflow::LoadSavedModel(tensorflow::SessionOptions(), + tensorflow::RunOptions(), + graph_file_name, + {tensorflow::kSavedModelTagServe}, + fBundle); + std::cout << "tf_graph loaded SavedModelBundle with status: " << status.ToString() << std::endl; + if (!status.ok()) return; + + auto sig_map = fBundle->meta_graph_def.signature_def(); + std::string sig_def = "serving_default"; + bool has_default_key = false; + std::vector sig_map_keys; + for (auto const& p : sig_map) { + if (p.first == sig_def) has_default_key = true; + sig_map_keys.push_back(p.first); + } + auto model_def = sig_map.at((has_default_key) ? sig_def : sig_map_keys.back()); + + // ... Get the input names + if (inputs.empty()){ + std::cout << "tf_graph using all inputs:" << std::endl; + for (auto const& p : model_def.inputs()) { + fInputNames.push_back(p.second.name()); + std::cout << "tf_graph InputName: " << fInputNames.back() << std::endl; + std::cout << " key: " << p.first << " value: " << p.second.name() << std::endl; + } + } + else{ + std::cout << "tf_graph using selected inputs:" << std::endl; + for (const auto& s : inputs) { + for (auto const& p : model_def.inputs()) { + if (p.first == s) { + fInputNames.push_back(p.second.name()); + std::cout << " key: " << p.first << " value: " << p.second.name() << std::endl; + } + } + } + } + + // ... Get the output names + // .. get all outputs if no specific name provided + if (outputs.empty()) { + std::cout << "tf_graph using all outputs:" << std::endl; + for (auto const& p : model_def.outputs()) { + fOutputNames.push_back(p.second.name()); + std::cout << " key: " << p.first << " value: " << p.second.name() << std::endl; + } + } + // .. or use only the outputs whose keys are specified + else { + std::cout << "tf_graph using selected outputs:" << std::endl; + for (const auto& s : outputs) { + for (auto const& p : model_def.outputs()) { + if (p.first == s) { + fOutputNames.push_back(p.second.name()); + // std::cout << " key: " << p.first << " value: " << p.second.name() << std::endl; + } + } + } + } + if (fOutputNames.empty()) { + std::cout << "tf_graph did not find outputs in SaveModelBundle." << std::endl; + return; + } + } + else { + + tensorflow::GraphDef graph_def; + status = tensorflow::ReadBinaryProto(tensorflow::Env::Default(), graph_file_name, &graph_def); + std::cout << "tf_graph loaded ProtoBuf graph with status: " << status.ToString() << std::endl; + if (!status.ok()) return; + + size_t ng = graph_def.node().size(); + std::cout << "node size: " << ng << std::endl; + std::cout << "Graph input names: " << std::endl; + for (int i = 0; i < n_inputs; ++i) { + fInputNames.push_back(graph_def.node()[i].name()); + std::cout << graph_def.node()[i].name() << " - "; + } + + // // last node as output if no specific name provided + // if (outputs.empty()) { + // for (int i = n_outputs; i > 0; --i) { + // fOutputNames.push_back(graph_def.node()[ng - i].name()); + // } + // } + if (outputs.empty()){ + for (int i = n_inputs; i < graph_def.node().size(); i++){ + std::cout << " key: " << graph_def.node()[i].name() << std::endl; + fOutputNames.push_back(graph_def.node()[i].name()); + } + } + else // or last nodes with names containing provided strings + { + std::string last, current, basename, name; + for (size_t n = 0; n < ng; ++n) { + // name = graph_def.node()[n].name(); + name = graph_def.node(n).name(); + // std::cout << name << std::endl; + auto pos = name.find("/"); + if (pos != std::string::npos) { basename = name.substr(0, pos); } + else { + continue; + } + + bool found = false; + for (const auto& s : outputs) { + if (name.find(s) != std::string::npos) { + found = true; + break; + } + } + if (found) { + if (!last.empty() ) { fOutputNames.push_back(last); } //&& (basename != current) (estava no ifs) + current = basename; + last = name; + std::cout << "base: " << basename << " last: " << last << std::endl; + } + } + if (!last.empty()) { fOutputNames.push_back(last); } + } + if (fOutputNames.empty()) { + std::cout << "Output nodes not found in the graph." << std::endl; + return; + } + // else{ + // for (auto &k : fOutputNames) std::cout << "Output: " << k << std::endl; + // } + status = fSession->Create(graph_def); + if (!status.ok()) { + std::cout << status.ToString() << std::endl; + return; + } + } + + success = true; // ok, graph loaded from the file +} + +tf::Graph::~Graph() +{ + fSession->Close().IgnoreError(); + delete fSession; + if (fUseBundle) { delete fBundle; } +} +// ------------------------------------------------------------------- + +std::vector> tf::Graph::run(const std::vector>& x) +{ + if (x.empty() || x.front().empty()) { return std::vector>(); } + + long long int rows = x.size(), cols = x.front().size(); + + std::vector _x; + _x.push_back( + tensorflow::Tensor(tensorflow::DT_FLOAT, tensorflow::TensorShape({1, rows, cols, 1}))); + auto input_map = _x[0].tensor(); + + for (long long int r = 0; r < rows; ++r) { + const auto& row = x[r]; + for (long long int c = 0; c < cols; ++c) { + input_map(0, r, c, 0) = row[c]; + } + } + + auto result = run(_x); + if (!result.empty()) { return result.front(); } + else { + return std::vector>(); + } +} +// ------------------------------------------------------------------- + +std::vector>> tf::Graph::run( + const std::vector>>>& x, + long long int samples) +{ + if ((samples == 0) || x.empty() || x.front().empty() || x.front().front().empty() || + x.front().front().front().empty()) + return std::vector>>(); + + if ((samples == -1) || (samples > (long long int)x.size())) { samples = x.size(); } + + long long int rows = x.front().size(), cols = x.front().front().size(), + depth = x.front().front().front().size(); + + std::vector _x; + + // Single-input network + if (n_inputs == 1) { + _x.push_back(tensorflow::Tensor(tensorflow::DT_FLOAT, + tensorflow::TensorShape({samples, rows, cols, depth}))); + auto input_map = _x[0].tensor(); + for (long long int s = 0; s < samples; ++s) { + const auto& sample = x[s]; + for (long long int r = 0; r < rows; ++r) { + const auto& row = sample[r]; + for (long long int c = 0; c < cols; ++c) { + const auto& col = row[c]; + for (long long int d = 0; d < depth; ++d) { + input_map(s, r, c, d) = col[d]; + } + } + } + } + } + // Multi-input network + else { + for (int i = 0; i < depth; ++i) { + _x.push_back(tensorflow::Tensor(tensorflow::DT_FLOAT, + tensorflow::TensorShape({samples, rows, cols, 1}))); + } + + for (int view = 0; view < depth; ++view) { + auto input_map = _x[view].tensor(); + for (long long int s = 0; s < samples; ++s) { + const auto& sample = x[s]; + for (long long int r = 0; r < rows; ++r) { + const auto& row = sample[r]; + for (long long int c = 0; c < cols; ++c) { + const auto& col = row[c]; + long long int d = view; + input_map(s, r, c, 0) = col[d]; + } + } + } + } + } + + return run(_x); +} + +// ------------------------------------------------------------------- + +std::vector>> tf::Graph::run( + const std::vector& x) +{ + std::vector> inputs; + for (int i = 0; i < n_inputs; ++i) { + inputs.push_back({fInputNames[i], x[i]}); + } + + std::vector outputs; + std::vector outputNames; + auto status = (fUseBundle) ? + fBundle->GetSession()->Run(inputs, fOutputNames, outputNames, &outputs) : + fSession->Run(inputs, fOutputNames, outputNames, &outputs); + + if (status.ok()) { + size_t samples = 0; + + for (size_t o = 0; o < outputs.size(); ++o) { + if (o == 0) { samples = outputs[o].dim_size(0); } + else if ((int)samples != outputs[o].dim_size(0)) { + throw std::string("TF outputs size inconsistent."); + } + } + + std::vector>> result; + result.resize(samples, std::vector>(outputs.size())); + + for (size_t s = 0; s < samples; ++s) { + for (size_t o = 0; o < outputs.size(); ++o) { + size_t n = outputs[o].dim_size(1); + auto output_map = outputs[o].tensor(); + + result[s][o].resize(outputs[o].dim_size(1)); + + std::vector& vs = result[s][o]; + for (size_t i = 0; i < n; ++i) { + vs[i] = output_map(s, i); + } + } + } + + return result; + } + else { + std::cout << status.ToString() << std::endl; + return std::vector>>(); + } +} +// ------------------------------------------------------------------- + +std::vector> tf::Graph::runx(const std::vector& x) +{ + std::vector> inputs; + for (int i = 0; i < n_inputs; ++i) { + inputs.push_back({fInputNames[i], x[i]}); + } + + std::vector outputs; + std::vector outputNames; + auto status = (fUseBundle) ? + fBundle->GetSession()->Run(inputs, fOutputNames, outputNames, &outputs) : + fSession->Run(inputs, fOutputNames, outputNames, &outputs); + + if (status.ok()) { + size_t samples = 0, nouts = 0; + + for (size_t o = 0; o < outputs.size(); ++o) { + if (o == 0) { samples = outputs[o].dim_size(0); } + else if ((int)samples != outputs[o].dim_size(0)) { + throw std::string("TF outputs size inconsistent."); + } + nouts += outputs[o].dim_size(1); + } + + std::vector> result; + result.resize(samples, std::vector(nouts)); + + size_t idx0 = 0; + for (size_t o = 0; o < outputs.size(); ++o) { + auto output_map = outputs[o].tensor(); + + size_t n = outputs[o].dim_size(1); + for (size_t s = 0; s < samples; ++s) { + std::vector& vs = result[s]; + for (size_t i = 0; i < n; ++i) { + vs[idx0 + i] = output_map(s, i); + } + } + idx0 += n; + } + + return result; + } + else { + std::cout << status.ToString() << std::endl; + return std::vector>(); + } +} +// ------------------------------------------------------------------- + +std::vector> tf::Graph::runae(const std::vector& x) +{ + std::vector> inputs; + for (int i = 0; i < n_inputs; ++i) { + inputs.push_back({fInputNames[i], x[i]}); + } + + std::vector outputs; + std::vector outputNames; + auto status = (fUseBundle) ? + fBundle->GetSession()->Run(inputs, fOutputNames, outputNames, &outputs) : + fSession->Run(inputs, fOutputNames, outputNames, &outputs); + + if (status.ok()) { + size_t samples = 0, npoints = 0; + + if (outputs.size() > 1) { throw std::string("TF runae: detected more than one output."); } + + samples = outputs[0].dim_size(0); + npoints = outputs[0].dim_size(1); + + std::vector> result; + result.resize(samples, std::vector(npoints)); + + auto output_map = outputs[0].tensor(); + + for (size_t s = 0; s < samples; ++s) { + std::vector& vs = result[s]; + for (size_t i = 0; i < npoints; ++i) { + vs[i] = output_map(s, i, 0); + } + } + + return result; + } + else { + std::cout << status.ToString() << std::endl; + return std::vector>(); + } +} diff --git a/icaruscode/ICARUSCVN/tf/tf_graph.h b/icaruscode/ICARUSCVN/tf/tf_graph.h new file mode 100644 index 000000000..0dd52ac7f --- /dev/null +++ b/icaruscode/ICARUSCVN/tf/tf_graph.h @@ -0,0 +1,80 @@ +//////////////////////////////////////////////////////////////////////////////////////////////////// +//// Class: Graph +//// Authors: R.Sulej (Robert.Sulej@cern.ch), from DUNE, FNAL/NCBJ, Sept. 2017 +//// P.Plonski, from DUNE, WUT, Sept. 2017 +//// T.Cai (tejinc@yorku.ca) from DUNE, YorkU, March 2022 +//// B.N.Nayak (nayakb@uci.edu) from DUNE, BNL, November 2022 +//// +//// Iterface to run Tensorflow graph saved to a file. First attempts, almost functional. +//// +//////////////////////////////////////////////////////////////////////////////////////////////////// + +#ifndef Graph_h +#define Graph_h + +#include +#include +#include + +namespace tensorflow { + class Session; + class Tensor; + struct SavedModelBundle; +} + +namespace tf { + + class Graph { + public: + static std::unique_ptr create(const char* graph_file_name, + const std::vector& inputs = {}, + const std::vector& outputs = {}, + bool use_bundle = false, + int ninputs = 1, + int noutputs = 1) + { + bool success; + std::unique_ptr ptr( + new Graph(graph_file_name, inputs, outputs, success, use_bundle, ninputs, noutputs)); + if (success) { return ptr; } + else { + return nullptr; + } + } + + ~Graph(); + + std::vector> run(const std::vector>& x); + + // process vector of 3D inputs, return vector of 1D outputs; use all inputs + // if samples = -1, or only the specified number of first samples; + // can deal with multiple inputs + std::vector>> run( + const std::vector>>>& x, + long long int samples = -1); + std::vector>> run(const std::vector& x); + std::vector> runx(const std::vector& x); + std::vector> runae(const std::vector& x); + + private: + int n_inputs; + int n_outputs; + /// Not-throwing constructor. + Graph(const char* graph_file_name, + const std::vector& inputs, + const std::vector& outputs, + bool& success, + bool use_bundle = false, + int ninputs = 1, + int noutputs = 1); + + tensorflow::Session* fSession; + bool fUseBundle; + tensorflow::SavedModelBundle* fBundle; + std::vector fInputNames; + std::vector fOutputNames; + }; + +} // namespace tf + +#endif From 8ac51f164a8c2e1fea0d9cf475200e80a873657b Mon Sep 17 00:00:00 2001 From: Felipe Wieler Date: Tue, 18 Nov 2025 08:39:51 -0600 Subject: [PATCH 2/2] Adding CVN code --- icaruscode/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/icaruscode/CMakeLists.txt b/icaruscode/CMakeLists.txt index 606f87f42..b857124fb 100644 --- a/icaruscode/CMakeLists.txt +++ b/icaruscode/CMakeLists.txt @@ -15,4 +15,5 @@ add_subdirectory(RecoUtils) add_subdirectory(TPC) add_subdirectory(Utilities) add_subdirectory(Overlays) +add_subdirectory(ICARUSCVN) #add_subdirectory(Supera)