From ae3374a03f98358c51ac3d76a4cb4afab0806ce3 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Fri, 18 Oct 2024 11:03:37 -0500 Subject: [PATCH 1/9] add module exporting visibility map and plotting/converting utils --- CMakeLists.txt | 3 + duneopdet/PhotonPropagation/CMakeLists.txt | 105 ++-- .../PhotonVisibilityExport.fcl | 25 + .../PhotonVisibilityExport_module.cc | 488 ++++++++++++++++++ .../VisibilityMapTools/ttree_to_th3.C | 144 ++++++ 5 files changed, 717 insertions(+), 48 deletions(-) create mode 100644 duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl create mode 100644 duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C diff --git a/CMakeLists.txt b/CMakeLists.txt index a11dcd59..1131ed0d 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,9 @@ find_ups_root() find_ups_product( cetbuildtools ) find_ups_product( larevt ) find_ups_product( larsim ) +find_ups_product( eigen ) +find_ups_product( tensorflow ) +find_ups_product( larsimdnn ) find_ups_product( larcore ) find_ups_product( lardata ) find_ups_product( lardataalg ) diff --git a/duneopdet/PhotonPropagation/CMakeLists.txt b/duneopdet/PhotonPropagation/CMakeLists.txt index 9d063eeb..f58a974a 100644 --- a/duneopdet/PhotonPropagation/CMakeLists.txt +++ b/duneopdet/PhotonPropagation/CMakeLists.txt @@ -1,54 +1,63 @@ art_make( - SERVICE_LIBRARIES larsim_PhotonPropagation - larsim::Simulation nug4::ParticleNavigation - lardataobj::Simulation - larevt::Filters - lardataobj::RawData - larcorealg::Geometry - larcore::Geometry_Geometry_service - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art_root_io::tfile_support - ROOT::Core - art_root_io::TFileService_service - art::Framework_Services_Optional_RandomNumberGenerator_service - art::Persistency_Common - art::Persistency_Provenance - art::Utilities - canvas::canvas - messagefacility::MF_MessageLogger - fhiclcpp::fhiclcpp - CLHEP - cetlib::cetlib cetlib_except - ROOT_BASIC_LIB_LIST - ROOT_EG + SERVICE_LIBRARIES + larsim_PhotonPropagation + larsim::Simulation nug4::ParticleNavigation + lardataobj::Simulation + larevt::Filters + lardataobj::RawData + larcorealg::Geometry + larcore::Geometry_Geometry_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support + ROOT::Core + art_root_io::TFileService_service + art::Framework_Services_Optional_RandomNumberGenerator_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + CLHEP + cetlib::cetlib cetlib_except + ROOT_BASIC_LIB_LIST + ROOT_EG + MODULE_LIBRARIES - ROOT_BASIC_LIB_LIST - larsim::LegacyLArG4 - lardataobj::Simulation - duneopdet::PhotonPropagation_PhotonVisibilityServiceS2_service - larsim::Simulation - nug4::ParticleNavigation lardataobj_Simulation - larcorealg::Geometry - larcore::Geometry_Geometry_service - art::Framework_Core - art::Framework_Principal - art::Framework_Services_Registry - art_root_io::tfile_support - ROOT::Core - art_root_io::TFileService_service - art::Framework_Services_Optional_RandomNumberGenerator_service - nurandom::RandomUtils_NuRandomService_service - art::Persistency_Common - art::Persistency_Provenance - art::Utilities - canvas::canvas - messagefacility::MF_MessageLogger - fhiclcpp::fhiclcpp - cetlib::cetlib cetlib_except - CLHEP + + ROOT_BASIC_LIB_LIST + larsim::LegacyLArG4 + lardataobj::Simulation + duneopdet::PhotonPropagation_PhotonVisibilityServiceS2_service + larsim::PhotonPropagation_PhotonVisibilityService_service + larsim::Simulation + larsim::PhotonPropagation + TensorFlow::cc + TensorFlow::framework + larsimdnn::TFLoaderTool + + nug4::ParticleNavigation lardataobj_Simulation + larcorealg::Geometry + larcore::Geometry_Geometry_service + art::Framework_Core + art::Framework_Principal + art::Framework_Services_Registry + art_root_io::tfile_support + ROOT::Core + art_root_io::TFileService_service + art::Framework_Services_Optional_RandomNumberGenerator_service + nurandom::RandomUtils_NuRandomService_service + art::Persistency_Common + art::Persistency_Provenance + art::Utilities + canvas::canvas + messagefacility::MF_MessageLogger + fhiclcpp::fhiclcpp + cetlib::cetlib cetlib_except + CLHEP ) install_headers() diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl new file mode 100644 index 00000000..450c714d --- /dev/null +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl @@ -0,0 +1,25 @@ +#include "PDFastSim_dune.fcl" + +BEGIN_PROLOG + +standard_photovis: +{ + module_type: "PhotonVisibilityExport" + photovis_label: "photovis" + photovis_map: "photovis" + voxel_dx: "10.00" # voxel x mesh step (in cm) + voxel_dy: "10.00" # voxel y mesh step (in cm) + voxel_dz: "10.00" # voxel z mesh step (in cm) + vis_model: "semianalytical" # visibility model to be used + # (pick one between "semianalytical", "compgraph") + do_refl: @local::dunefd_pdfastsim_par_ar_refl.DoReflectedLight + do_include_anode_refl: "false" + vuvhitspars: @local::dunefd_pdfastsim_par_ar.VUVHits + vishitspars: @local::dunefd_pdfastsim_par_ar_refl.VISHits + tfloaderpars: @local::dunevd_pdfastsim_ann_ar.TFLoaderTool + + do_include_buffer: false +} + +END_PROLOG + diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc new file mode 100644 index 00000000..ffaf5de3 --- /dev/null +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc @@ -0,0 +1,488 @@ +/** + * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) + * @file : PhotonVisibilityExport_module.cc + * @created : Thursday Feb 10, 2022 05:19:33 CST + */ + +#ifndef PHOTONVISIBILITYEXPORT_MODULE_CC + +#define PHOTONVISIBILITYEXPORT_MODULE_CC + +// ROOT includes +#include "TH1D.h" +#include "TEfficiency.h" +#include "TFile.h" +#include "TTree.h" +#include "TVectorF.h" +#include "TRandom3.h" +#include "TParameter.h" + +// C++ includes +#include +#include +#include +#include +#include +#include +#include "math.h" +#include + +// LArSoft includes +#include "larcore/Geometry/Geometry.h" +#include "larsim/PhotonPropagation/SemiAnalyticalModel.h" +#include "larsimdnn/PhotonPropagation/TFLoaderTools/TFLoader.h" +#include "larsim/PhotonPropagation/PhotonVisibilityService.h" +#include "larcorealg/CoreUtils/counter.h" + +// ART includes. +#include "art/Framework/Core/EDAnalyzer.h" +#include "art/Framework/Principal/Event.h" +#include "fhiclcpp/ParameterSet.h" +#include "art/Framework/Principal/Handle.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 "messagefacility/MessageLogger/MessageLogger.h" +#include "art/Framework/Core/ModuleMacros.h" +#include "art/Utilities/make_tool.h" +#include "canvas/Persistency/Common/FindManyP.h" + +namespace opdet { + + class PhotonVisibilityExport : public art::EDAnalyzer + { + public: + enum EVisModel {kSemiAnalytical = 0, kCompGraph = 1}; + + PhotonVisibilityExport(const fhicl::ParameterSet&); + virtual ~PhotonVisibilityExport() {}; + + void beginJob(); + void endJob() {} + void analyze (const art::Event&); + + private: + std::unique_ptr fVisibilityModel; + std::unique_ptr fTFGenerator; + size_t fNOpDets; + std::vector fOpDetCenter; + + EVisModel kVisModel; + + bool fDoReflectedLight = {}; + bool fIncludeAnodeReflections = {}; + bool fIncludeBuffer = {}; + bool fUseXeAbsorption = {}; + + double fVoxelSizeX = {}; + double fVoxelSizeY = {}; + double fVoxelSizeZ = {}; + + fhicl::ParameterSet fVUVHitsParams; + fhicl::ParameterSet fVISHitsParams; + fhicl::ParameterSet fTFLoaderPars; + + bool fIsDone = false; + + std::array fCryostatMin = {}; + std::array fCryostatMax = {}; + std::array fTPCMin = {}; + std::array fTPCMax = {}; + + TH1D* fhGrid[3] = {}; + + void ExportTPCMap(); + void ExportOpDetMap(); + void ExportVoxelGrid(); + void ExportVisibility(); + + }; + +} // close opdet namespace + +namespace opdet { + DEFINE_ART_MODULE(PhotonVisibilityExport) +} + + +#endif /* end of include guard VISMAPDUMP_MODULE_CC */ + +namespace opdet { + PhotonVisibilityExport::PhotonVisibilityExport(fhicl::ParameterSet const& pset) : + art::EDAnalyzer(pset) + { + //* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * + // Read inputs from fcl file + fVoxelSizeX = pset.get("voxel_dx", 10.0); + fVoxelSizeY = pset.get("voxel_dy", 10.0); + fVoxelSizeZ = pset.get("voxel_dz", 10.0); + + fDoReflectedLight = pset.get("do_refl", false); + fIncludeAnodeReflections = pset.get("do_include_anode_refl", false); + fIncludeBuffer = pset.get("do_include_buffer", false); + fUseXeAbsorption = pset.get("do_include_xe_absorption", false); + + fVUVHitsParams = pset.get("vuvhitspars"); + fVISHitsParams = pset.get("vishitspars"); + + fTFLoaderPars = pset.get("tfloaderpars"); + + TString vis_model_str = pset.get("vis_model"); + + if (vis_model_str.Contains("compgraph")) { + kVisModel = kCompGraph; + } + else if (vis_model_str.Contains("semianalytical")) { + kVisModel = kSemiAnalytical; + } + } + + void PhotonVisibilityExport::beginJob() { + // create the photo-detector visibility model + if (kVisModel == kSemiAnalytical) { + printf("Creating Semi-analytical visibility model\n"); + fVisibilityModel = std::make_unique( + fVUVHitsParams, fVISHitsParams, + fDoReflectedLight, fIncludeAnodeReflections, fUseXeAbsorption + ); + } + else if (kVisModel == kCompGraph) { + fTFGenerator = art::make_tool(fTFLoaderPars); + fTFGenerator->Initialization(); + } + + return; + } + + void PhotonVisibilityExport::analyze(const art::Event&) { + ExportOpDetMap(); + + ExportTPCMap(); + + ExportVisibility(); + + return; + } + + void PhotonVisibilityExport::ExportVoxelGrid() { + art::ServiceHandle tfs; + + const double voxelDim[3] = {fVoxelSizeX, fVoxelSizeY, fVoxelSizeZ}; + // define mesh points with an histogram helper to align tpcs and buffer + // sampling coordinates + for (int i=0; i<3; i++) { + double xc = 0.5*(fCryostatMax[i] + fCryostatMin[i]); + int nbin = ceil( (fCryostatMax[i]-fCryostatMin[i]) / voxelDim[i] ); + double xmin = xc - nbin*0.5*voxelDim[i]; + double xmax = xc + nbin*0.5*voxelDim[i]; + fhGrid[i] = tfs->make(Form("hgrid%i", i), Form("mesh points axis %i", i), + nbin, xmin, xmax); + } + } + + void PhotonVisibilityExport::ExportTPCMap() { + // retrieve geometry + const auto geom = art::ServiceHandle(); + + art::ServiceHandle< art::TFileService > tfs; + + float tpcH = 0; + float tpcW = 0; + float tpcL = 0; + float tpcPos[3]; + //const double voxelDim[3] = {fVoxelSizeX, fVoxelSizeY, fVoxelSizeZ}; + + TTree* tTPC = tfs->make("tpcMap", "tpcMap"); + tTPC->Branch("tpcH", &tpcH); + tTPC->Branch("tpcL", &tpcL); + tTPC->Branch("tpcW", &tpcW); + tTPC->Branch("tpcPos", &tpcPos, "tpcPos[3]/F"); + + // get cryostat dimensions + for (geo::CryostatGeo const& cryo : geom->Iterate()) { + printf("cryostat: %u [%g, %g, %g] - [%g, %g, %g]\n", cryo.ID().getIndex(), + cryo.MinX(), cryo.MinY(), cryo.MinZ(), cryo.MaxX(), cryo.MaxY(), cryo.MaxZ()); + fCryostatMin[0] = cryo.MinX(); fCryostatMax[0] = cryo.MaxX(); + fCryostatMin[1] = cryo.MinY(); fCryostatMax[1] = cryo.MaxY(); + fCryostatMin[2] = cryo.MinZ(); fCryostatMax[2] = cryo.MaxZ(); + } + + + // loop over all TPCs to get the active volume dimensions + for (geo::TPCGeo const& tpc : geom->Iterate()) { + auto point_center = tpc.GetCenter(); + //point_center.GetCoordinates( tpcPos ); + double hlfW = tpc.ActiveHalfWidth (); tpcW = 2*hlfW; + double hlfH = tpc.ActiveHalfHeight(); tpcH = 2*hlfH; + double hlfL = tpc.ActiveHalfLength(); tpcL = 2*hlfL; + + const double hlfDim[3] = {hlfW, hlfH, hlfL}; + + double x_min[3] = {0}; double x_max[3] = {0}; + point_center.GetCoordinates( x_min ); + point_center.GetCoordinates( x_max ); + for (int idim =0; idim<3; idim++) { + x_min[idim] -= hlfDim[idim]; + x_max[idim] += hlfDim[idim]; + + if (fTPCMin[idim] > x_min[idim]) fTPCMin[idim] = x_min[idim]; + if (fTPCMax[idim] < x_max[idim]) fTPCMax[idim] = x_max[idim]; + } + tTPC->Fill(); + + printf("TPC dimensions: hlfW %.2f - hlfH %.2f hlfL %.2f, ", + hlfW, hlfH, hlfL); + printf("TPC center: %.2f - %.2f - %.2f\n", + point_center.x(), point_center.y(), point_center.z()); + } + + printf("Cryostat min-max: [%g, %g, %g] - [%g, %g, %g]\n", + fCryostatMin[0], fCryostatMin[1], fCryostatMin[2], + fCryostatMax[0], fCryostatMax[1], fCryostatMax[2]); + printf("TPC min-max: [%g, %g, %g] - [%g, %g, %g]\n", + fTPCMin[0], fTPCMin[1], fTPCMin[2], + fTPCMax[0], fTPCMax[1], fTPCMax[2]); + + TString descriptor = ""; + Double_t tmp[3]; + + TTree* tDimensions = tfs->make("tDimensions", "TPC and cryostat dimensions"); + tDimensions->Branch("descriptor", &descriptor); + tDimensions->Branch("dimension", &tmp, "dimension[3]/D"); + + descriptor = "tpc_min"; + for (size_t i = 0; i < 3; i++) tmp[i] = fTPCMin[i]; + tDimensions->Fill(); + + descriptor = "tpc_max"; + for (size_t i = 0; i < 3; i++) tmp[i] = fTPCMax[i]; + tDimensions->Fill(); + + descriptor = "cryostat_min"; + for (size_t i = 0; i < 3; i++) tmp[i] = fCryostatMin[i]; + tDimensions->Fill(); + + descriptor = "cryostat_max"; + for (size_t i = 0; i < 3; i++) tmp[i] = fCryostatMax[i]; + tDimensions->Fill(); + + ExportVoxelGrid(); + + return; + } + + void PhotonVisibilityExport::ExportOpDetMap() { + // retrieve geometry + const auto geom = art::ServiceHandle(); + + art::ServiceHandle< art::TFileService > tfs; + + float opDetH = 0; + float opDetW = 0; + float opDetL = 0; + float opDetPos[3]; + std::vector opChannel; + + TTree* tOpDet = tfs->make("opDetMap", "opDetMap"); + tOpDet->Branch("opDetH", &opDetH); + tOpDet->Branch("opDetL", &opDetL); + tOpDet->Branch("opDetW", &opDetW); + tOpDet->Branch("opDetPos", &opDetPos, "opDetPos[3]/F"); + tOpDet->Branch("opDetCh", &opChannel); + + // store info from Geometry service + fNOpDets = geom->NOpDets(); + for (size_t i : util::counter(fNOpDets)) { + opChannel.clear(); + geo::OpDetGeo const& opDet = geom->OpDetGeoFromOpDet(i); + auto center = opDet.GetCenter(); + center.GetCoordinates( opDetPos ); + opDetH = opDet.Height(); + opDetW = opDet.Width(); + opDetL = opDet.Length(); + + size_t n_ch = geom->NOpHardwareChannels(i); + opChannel.resize(n_ch, 0); + for (size_t ich = 0; ich < n_ch; ich++) { + opChannel.at(ich) = geom->OpChannel(i, ich); + } + + tOpDet->Fill(); + } + return; + } + + + void PhotonVisibilityExport::ExportVisibility() { + const auto photonVisService = art::ServiceHandle(); + // open file + art::ServiceHandle< art::TFileService > tfs; + + Double_t point_[3]; + double total_visDirect = 0; + double total_visReflct = 0; + double* opDet_visDirect = new double[fNOpDets]; + double* opDet_visReflct = new double[fNOpDets]; + Double_t pointBuff_[3]; + double total_visDirectBuff = 0; + double total_visReflctBuff = 0; + double* opDet_visDirectBuff = new double[fNOpDets]; + double* opDet_visReflctBuff = new double[fNOpDets]; + + // open tree + TTree* tMap = tfs->make("photoVisMap", "photoVisMap"); + tMap->SetNameTitle("photoVisMap", "photoVisMap"); + tMap->Branch("coords", &point_, "coords[3]/D"); + tMap->Branch("opDet_visDirect", opDet_visDirect, Form("opDet_visDirect[%ld]/D", fNOpDets)); + tMap->Branch("opDet_visReflct", opDet_visReflct, Form("opDet_visReflct[%ld]/D", fNOpDets)); + tMap->Branch("total_visDirect", &total_visDirect, "total_visDirect/D"); + tMap->Branch("total_visReflct", &total_visReflct, "total_visReflct/D"); + + TTree* tMapBuffer = tfs->make("photoVisMapBuffer", "photoVisMapBuffer"); + tMapBuffer->SetNameTitle("photoVisMapBuffer", "photoVisMapBuffer"); + tMapBuffer->Branch("coords", &pointBuff_, "coords[3]/D"); + tMapBuffer->Branch("opDet_visDirectBuff", opDet_visDirectBuff, Form("opDet_visDirectBuff[%ld]/D", fNOpDets)); + tMapBuffer->Branch("opDet_visReflctBuff", opDet_visReflctBuff, Form("opDet_visReflctBuff[%ld]/D", fNOpDets)); + tMapBuffer->Branch("total_visDirectBuff", &total_visDirectBuff, "total_visDirectBuff/D"); + tMapBuffer->Branch("total_visReflctBuff", &total_visReflctBuff, "total_visReflctBuff/D"); + + std::vector opdetvis_dir; + std::vector opdetvis_rfl; + + bool tpc_range_x = false, tpc_range_y = false, tpc_range_z = false; + const double voxelDim[3] = {fVoxelSizeX, fVoxelSizeY, fVoxelSizeZ}; + + // here we can loop over the points of the pre-defined grid + double x_= 0, y_= 0, z_= 0; + printf("starting loop: fHGrid histos: [%p, %p, %p]\n", + static_cast(fhGrid[0]), static_cast(fhGrid[1]), static_cast(fhGrid[2])); + + for (int ix=1; ix<=fhGrid[0]->GetNbinsX(); ix++) { + x_ = fhGrid[0]->GetBinCenter(ix); + if (x_> fTPCMin[0] && x_< fTPCMax[0]) tpc_range_x = true; + else tpc_range_x = false; + + for (int iy=1; iy<=fhGrid[1]->GetNbinsX(); iy++) { + y_= fhGrid[1]->GetBinCenter(iy); + if (y_> fTPCMin[1] && y_< fTPCMax[1]) tpc_range_y = true; + else tpc_range_y = false; + + for (int iz=1; iz<=fhGrid[2]->GetNbinsX(); iz++) { + z_= fhGrid[2]->GetBinCenter(iz); + if (z_> fTPCMin[2] && z_< fTPCMax[2]) tpc_range_z = true; + else tpc_range_z = false; + + opdetvis_dir.clear(); + opdetvis_rfl.clear(); + + if ( (tpc_range_x && tpc_range_y && tpc_range_z) == false) { + if (fIncludeBuffer) { + total_visDirectBuff = 0; + total_visReflctBuff = 0; + for (size_t i = 0; i < fNOpDets; i++) { + opDet_visDirectBuff[i] = 0.; + opDet_visReflctBuff[i] = 0.; + } + + const geo::Point_t point( x_, y_, z_); + auto mapped_vis = photonVisService->GetAllVisibilities(point); + size_t iopdet = 0; + for (auto &vis : mapped_vis) { + if (vis < 0.0) { + fprintf(stderr, "WARNIGN: negative buffer vis at (%.0f, %.0f, %.0f). visibility is %g\n", + x_, y_, z_, vis); + } + else if (vis > 1.0) { + fprintf(stderr, "WARNIGN: tpc vis at (%.0f, %.0f, %.0f) is greater than 1. visibility is %g\n", + x_, y_, z_, vis); + } + else { + opDet_visDirectBuff[iopdet] = vis; + total_visDirectBuff += vis; + } + iopdet++; + } + + point.GetCoordinates( pointBuff_ ); + tMapBuffer->Fill(); + } + } + else { + const geo::Point_t vpoint( x_, y_, z_); + total_visDirect = 0.; + total_visReflct = 0.; + for (size_t i = 0; i < fNOpDets; i++) { + opDet_visDirect[i] = 0.0; + opDet_visReflct[i] = 0.0; + } + + const double n_samplings = 5.0; + + for (int i = 0; i < n_samplings; i++) { + const geo::Vector_t delta( + gRandom->Uniform(-0.5*voxelDim[0], 0.5*voxelDim[0]), + gRandom->Uniform(-0.5*voxelDim[1], 0.5*voxelDim[1]), + gRandom->Uniform(-0.5*voxelDim[2], 0.5*voxelDim[2]) ); + + const geo::Point_t xspot = vpoint + delta; + + if (kVisModel == kSemiAnalytical ) { + fVisibilityModel->detectedDirectVisibilities (opdetvis_dir, xspot); + //fVisibilityModel->detectedReflectedVisibilities(opdetvis_rfl, xspot, fIncludeAnodeReflections); + } + else if (kVisModel == kCompGraph ) { + std::vector pos_tmp {xspot.x(), xspot.y(), xspot.z()}; + fTFGenerator->Predict( pos_tmp ); + opdetvis_dir = fTFGenerator->GetPrediction(); + } + + size_t iopdet = 0; + for (const auto &vis : opdetvis_dir) { + if (vis < 0.0) { + fprintf(stderr, "WARNIGN: negative tpc vis at (%.0f, %.0f, %.0f). visibility is %g\n", + x_, y_, z_, vis); + } + else { + opDet_visDirect[iopdet] += vis; + total_visDirect += vis; + } + iopdet++; + } + + iopdet = 0; + for (const auto &vis : opdetvis_rfl) { + opDet_visReflct[iopdet] += vis; + total_visReflct += vis; + iopdet++; + } + } + + total_visDirect = total_visDirect / n_samplings; + //total_visReflct = total_visReflct / n_samplings; + + for (size_t i = 0; i < fNOpDets; i++) { + opDet_visDirect[i] /= n_samplings; + opDet_visReflct[i] /= n_samplings; + //printf("[%ld] - %g\n", i, opDet_visDirect[i]); + } + + // Fill tree + vpoint.GetCoordinates( point_ ); + tMap->Fill(); + } + } + } + } + + fIsDone = true; + delete[] opDet_visReflct; + delete[] opDet_visDirect; + delete[] opDet_visDirectBuff; + delete[] opDet_visReflctBuff; + return; + } + + } + diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C new file mode 100644 index 00000000..2874400d --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C @@ -0,0 +1,144 @@ +/** + * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) + * @file : ttree_to_th3.C + * @created : 2024-10-18 08:13 + */ + +#include +#include +#include +#include +#include + +#include "RtypesCore.h" +#include "TFile.h" +#include "TTree.h" +#include "TTreeReader.h" +#include "TTreeReaderArray.h" +#include "TTreeReaderValue.h" +#include "TDirectoryFile.h" +#include "TSystem.h" +#include "TStyle.h" +#include "TCanvas.h" +#include "TClonesArray.h" +#include "TH1D.h" +#include "TH3F.h" +#include +#include +#include "TH2D.h" +#include "TGraphErrors.h" +#include "TF1.h" +#include "TRandom3.h" +#include "TROOT.h" + +void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_t vis_threshold = 1e-4) +{ + gStyle->SetTitleSize(0.06, "xyz"); + gStyle->SetLabelSize(0.06, "xyz"); + + TFile* file = new TFile(input_file); + TTree* tOpDet = nullptr; + TTree* tVisTPC = nullptr; + TTree* tVisBuf = nullptr; + TDirectoryFile* dunevis_dir = file->Get(vis_dir); + dunevis_dir->cd(); + + TH1D* h_tmp[3] = {0}; + h_tmp[0] = dunevis_dir->Get("hgrid0"); + h_tmp[1] = dunevis_dir->Get("hgrid1"); + h_tmp[2] = dunevis_dir->Get("hgrid2"); + + const int nbins[3] = {h_tmp[0]->GetNbinsX(), h_tmp[1]->GetNbinsX(), h_tmp[2]->GetNbinsX()}; + const double xmin[3] ={ + h_tmp[0]->GetXaxis()->GetXmin(), + h_tmp[1]->GetXaxis()->GetXmin(), + h_tmp[2]->GetXaxis()->GetXmin()}; + const double xmax[3] ={ + h_tmp[0]->GetXaxis()->GetXmax(), + h_tmp[1]->GetXaxis()->GetXmax(), + h_tmp[2]->GetXaxis()->GetXmax()}; + + tOpDet = dunevis_dir->Get("opDetMap"); + tVisTPC = dunevis_dir->Get("photoVisMap"); + tVisBuf = dunevis_dir->Get("photoVisMapBuffer"); + + TTreeReader reader(tVisTPC); + TTreeReaderArray _xyz(reader, "coords"); + TTreeReaderArray _opdet_visd(reader, "opDet_visDirect"); + //TTreeReaderArray _opdet_visr(reader, "opDet_visReflct"); + TTreeReaderValue _vis_d(reader, "total_visDirect"); + //TTreeReaderValue _vis_r(reader, "total_visReflct"); + + const Long64_t N_OPDET = tOpDet->GetEntries(); + + THnF* h3total_vis = new THnF("h3VisMap", "Visibility Map", 3, nbins, xmin, xmax); + std::vector h3opDet; + h3opDet.reserve(N_OPDET); + for (size_t iopdet = 0; iopdet < N_OPDET; iopdet++) { + TString name = Form("h3VisMap_opDet%ld", iopdet); + TString titl = Form("Visibility Map - OpDet %ld", iopdet); + h3opDet.emplace_back(new THnSparseF(name, titl, 3, nbins, xmin, xmax)); + } + + TH1D* h_opdet_vis_tpc = new TH1D("h_opdet_vis_tpc", "opdet visibility (TPC)", 100, -15, 1); + TH1D* h_opdet_vis_buf = new TH1D("h_opdet_vis_buf", "opdet visibility (buffer)", 100, -15, 1); + + while( reader.Next() ) { + double xyz[3] = {_xyz[0], _xyz[1], _xyz[2]}; + h3total_vis->Fill(xyz, *(_vis_d)); + for (Long64_t i = 0; i < N_OPDET; i++) { + h_opdet_vis_tpc->Fill( log10(_opdet_visd[i]) ); + if (_opdet_visd[i] > vis_threshold) { + h3opDet[i]->Fill(xyz, _opdet_visd[i]); + } + } + } + + if (tVisBuf) { + if (tVisBuf->GetEntries() > 0) { + TTreeReader readerBuff(tVisBuf); + TTreeReaderArray _buf_xyz(readerBuff, "coords"); + TTreeReaderArray _buf_opdet_visd(readerBuff, "opDet_visDirectBuff"); + //TTreeReaderArray _buf_opdet_visr(readerBuff, "opDet_visReflctBuff"); + TTreeReaderValue _buf_vis_d(readerBuff, "total_visDirectBuff"); + //TTreeReaderValue _buf_vis_r(readerBuff, "total_visReflctBuff"); + + while( readerBuff.Next() ) { + double xyz[3] = {_buf_xyz[0], _buf_xyz[1], _buf_xyz[2]}; + h3total_vis->Fill(xyz, *(_buf_vis_d)); + for (Long64_t i = 0; i < N_OPDET; i++) { + h_opdet_vis_buf->Fill( log10(_buf_opdet_visd[i]) ); + if ( _buf_opdet_visd[i] > vis_threshold ) { + h3opDet[i]->Fill(xyz, _buf_opdet_visd[i]); + } + } + } + } + } + + TString output_name = input_file; + output_name.Resize( output_name.Index(".root") ); + output_name += "_"+vis_dir+".root"; + + TFile* output_h3 = new TFile(output_name, "recreate"); + h3total_vis->Write(); + for (auto& h3 : h3opDet) { + h3->Write(); + } + output_h3->Close(); + + return; +} + +TH3F* rebin_visibility_map(const TH3F* h3, const int rbx, const int rby, const int rbz) +{ + TH3F* h3rb = (TH3F*)h3->Clone(); + h3rb = (TH3F*)h3rb->Rebin3D(rbx, rby, rbz, Form("%s_rb_%i_%i_%i", h3->GetName(), rbx, rby, rbz)); + + float scale_factor = 1.0 / (rbx * rby *rbz); + h3rb->Scale( scale_factor ); + + return h3rb; +} + + From 5d15e0917251b34adcaa08945adc4c38519478f6 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Tue, 22 Oct 2024 04:10:31 -0500 Subject: [PATCH 2/9] minor editing to projection function --- .../VisibilityMapTools/DUNEVisUtils.hpp | 149 ++++++++++++++++++ 1 file changed, 149 insertions(+) create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp new file mode 100644 index 00000000..d62b1183 --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp @@ -0,0 +1,149 @@ +/** + * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) + * @file : DUNEVisUtils.hpp + * @created : Friday Oct 18, 2024 17:19:58 CEST + */ + +#ifndef DUNEVISUTILS_HPP + +#define DUNEVISUTILS_HPP + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +inline TH3* rebin_visibility_map(const TH3* h3, const int rbx, const int rby, const int rbz) +{ + TH3* h3rb = (TH3*)h3rb->Rebin3D(rbx, rby, rbz, Form("%s_rb_%i_%i_%i", h3->GetName(), rbx, rby, rbz)); + + float scale_factor = 1.0 / (rbx * rby *rbz); + h3rb->Scale( scale_factor ); + return h3rb; +} + +inline THnSparse* rebin_visibility_map(const THnSparse* hn, const int rbx, const int rby, const int rbz) +{ + const Int_t rbidx[3] = {rbx, rby, rbz}; + THnSparse* hnrb = hn->Rebin(rbidx); + + float scale_factor = 1.0 / (rbx * rby * rbz); + hnrb->Scale( scale_factor ); + return hnrb; +} + +inline THn* rebin_visibility_map(const THn* hn, const int rbx, const int rby, const int rbz) { + const Int_t rbidx[3] = {rbx, rby, rbz}; + THn* hnrb = hn->Rebin(rbidx); + + float scale_factor = 1.0 / (rbx * rby * rbz); + hnrb->Scale( scale_factor ); + return hnrb; +} + +inline const int get_projection_axis(const int& axis0, const int& axis1) { + assert(axis0 != axis1); + assert(axis1 < 2); + assert(axis0 < 2); + + return (3-axis0-axis1); +} + +inline TH2D* draw_projection(const THnBase* hn, + const int& axis0, + const int& axis1, + const double& x0, + const double& x1) +{ + const int proj_axis = get_projection_axis(axis0, axis1); + const int ibin0 = hn->GetAxis(proj_axis)->FindBin(x0); + const int ibin1 = hn->GetAxis(proj_axis)->FindBin(x1); + const int proj_nbin = ibin1 - ibin0 + 1; + hn->GetAxis(proj_axis)->SetRange(ibin0, ibin1); + + TH2D* h_proj = hn->Projection(axis0, axis1); + h_proj->Scale( 1.0 / proj_nbin ); + + hn->GetAxis(proj_axis)->SetRange(); + return h_proj; +} + +inline void CanvasPartition(TCanvas *C,const Int_t Nx,const Int_t Ny, + Float_t lMargin, Float_t rMargin, + Float_t bMargin, Float_t tMargin) +{ + if (!C) return; + // Setup Pad layout: + Float_t vSpacing = 0.0; + Float_t vStep = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny; + Float_t hSpacing = 0.0; + Float_t hStep = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx; + Float_t vposd,vposu,vmard,vmaru,vfactor; + Float_t hposl,hposr,hmarl,hmarr,hfactor; + for (Int_t i=0;icd(0); + char name[16]; + sprintf(name,"pad_%i_%i",i,j); + TPad *pad = (TPad*) gROOT->FindObject(name); + if (pad) delete pad; + pad = new TPad(name,"",hposl,vposd,hposr,vposu); + pad->SetLeftMargin(hmarl); + pad->SetRightMargin(hmarr); + pad->SetBottomMargin(vmard); + pad->SetTopMargin(vmaru); + pad->SetFrameBorderMode(0); + pad->SetBorderMode(0); + pad->SetBorderSize(0); + pad->Draw(); + } + } +} + + +#endif /* end of include guard DUNEVISUTILS_HPP */ + From c499d6f3df06e50d475c51206ff53cb17dfa76bd Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Tue, 22 Oct 2024 08:07:48 -0500 Subject: [PATCH 3/9] Basic documentation --- .../PhotonVisibilityExport_module.cc | 59 ++++---- .../VisibilityMapTools/DUNEVisUtils.hpp | 128 ++++++------------ .../VisibilityMapTools/ttree_to_th3.C | 41 +++--- 3 files changed, 93 insertions(+), 135 deletions(-) diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc index ffaf5de3..f8e1100a 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc @@ -1,7 +1,7 @@ /** * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) * @file : PhotonVisibilityExport_module.cc - * @created : Thursday Feb 10, 2022 05:19:33 CST + * @created : 2024-10-22 07:50 */ #ifndef PHOTONVISIBILITYEXPORT_MODULE_CC @@ -51,6 +51,22 @@ namespace opdet { + /** + * @class PhotonVisibilityExport + * @brief Export the visibility map of the detector volume and of each individual OpDet + * + * This module produces a visibility map of the detector volume and of each individual + * optical detector, both inside the TPC volume (photoVisMap) and in the buffer region + * (photoVisMapBuffer). The TTree products can be converted into a 3D THn histogram with the + * macro in the VisibilityMapTools/ folder. + * In addition to the visibility tree, the module exports a tree containing the placement + * of the individual optical detector and a tree with the size and placementes of the + * detector's TPCs. Finally, the tDimensions tree can be used to retrieve the overall size + * of the detector. + * + * The voxel grid used for sampling the detector visibility can be recovered on the basis + * of the three histograms hgrid[0-2]. + */ class PhotonVisibilityExport : public art::EDAnalyzer { public: @@ -157,12 +173,14 @@ namespace opdet { } void PhotonVisibilityExport::analyze(const art::Event&) { - ExportOpDetMap(); - ExportTPCMap(); - - ExportVisibility(); + if (fIsDone == false) { + ExportOpDetMap(); + + ExportTPCMap(); + ExportVisibility(); + } return; } @@ -192,7 +210,6 @@ namespace opdet { float tpcW = 0; float tpcL = 0; float tpcPos[3]; - //const double voxelDim[3] = {fVoxelSizeX, fVoxelSizeY, fVoxelSizeZ}; TTree* tTPC = tfs->make("tpcMap", "tpcMap"); tTPC->Branch("tpcH", &tpcH); @@ -213,7 +230,6 @@ namespace opdet { // loop over all TPCs to get the active volume dimensions for (geo::TPCGeo const& tpc : geom->Iterate()) { auto point_center = tpc.GetCenter(); - //point_center.GetCoordinates( tpcPos ); double hlfW = tpc.ActiveHalfWidth (); tpcW = 2*hlfW; double hlfH = tpc.ActiveHalfHeight(); tpcH = 2*hlfH; double hlfL = tpc.ActiveHalfLength(); tpcL = 2*hlfL; @@ -390,18 +406,8 @@ namespace opdet { auto mapped_vis = photonVisService->GetAllVisibilities(point); size_t iopdet = 0; for (auto &vis : mapped_vis) { - if (vis < 0.0) { - fprintf(stderr, "WARNIGN: negative buffer vis at (%.0f, %.0f, %.0f). visibility is %g\n", - x_, y_, z_, vis); - } - else if (vis > 1.0) { - fprintf(stderr, "WARNIGN: tpc vis at (%.0f, %.0f, %.0f) is greater than 1. visibility is %g\n", - x_, y_, z_, vis); - } - else { - opDet_visDirectBuff[iopdet] = vis; - total_visDirectBuff += vis; - } + opDet_visDirectBuff[iopdet] = vis; + total_visDirectBuff += vis; iopdet++; } @@ -430,7 +436,7 @@ namespace opdet { if (kVisModel == kSemiAnalytical ) { fVisibilityModel->detectedDirectVisibilities (opdetvis_dir, xspot); - //fVisibilityModel->detectedReflectedVisibilities(opdetvis_rfl, xspot, fIncludeAnodeReflections); + fVisibilityModel->detectedReflectedVisibilities(opdetvis_rfl, xspot, fIncludeAnodeReflections); } else if (kVisModel == kCompGraph ) { std::vector pos_tmp {xspot.x(), xspot.y(), xspot.z()}; @@ -440,14 +446,8 @@ namespace opdet { size_t iopdet = 0; for (const auto &vis : opdetvis_dir) { - if (vis < 0.0) { - fprintf(stderr, "WARNIGN: negative tpc vis at (%.0f, %.0f, %.0f). visibility is %g\n", - x_, y_, z_, vis); - } - else { - opDet_visDirect[iopdet] += vis; - total_visDirect += vis; - } + opDet_visDirect[iopdet] += vis; + total_visDirect += vis; iopdet++; } @@ -460,12 +460,11 @@ namespace opdet { } total_visDirect = total_visDirect / n_samplings; - //total_visReflct = total_visReflct / n_samplings; + total_visReflct = total_visReflct / n_samplings; for (size_t i = 0; i < fNOpDets; i++) { opDet_visDirect[i] /= n_samplings; opDet_visReflct[i] /= n_samplings; - //printf("[%ld] - %g\n", i, opDet_visDirect[i]); } // Fill tree diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp index d62b1183..a93bd52f 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp @@ -19,34 +19,51 @@ #include #include -inline TH3* rebin_visibility_map(const TH3* h3, const int rbx, const int rby, const int rbz) -{ - TH3* h3rb = (TH3*)h3rb->Rebin3D(rbx, rby, rbz, Form("%s_rb_%i_%i_%i", h3->GetName(), rbx, rby, rbz)); - - float scale_factor = 1.0 / (rbx * rby *rbz); - h3rb->Scale( scale_factor ); - return h3rb; -} -inline THnSparse* rebin_visibility_map(const THnSparse* hn, const int rbx, const int rby, const int rbz) +/** + * @brief Rebin the visibility map averaging the bin content after rebin + * + * @param hn visibility map + * @param rb0 bin group for axis 0 + * @param rb1 bin group for axis 1 + * @param rb2 bin group for axis 2 + * @return rebinned visibility map + */ +inline THnSparse* rebin_visibility_map( + const THnSparse* hn, + const int rb0, + const int rb1, + const int rb2) { - const Int_t rbidx[3] = {rbx, rby, rbz}; + const Int_t rbidx[3] = {rb0, rb1, rb2}; THnSparse* hnrb = hn->Rebin(rbidx); - float scale_factor = 1.0 / (rbx * rby * rbz); + float scale_factor = 1.0 / (rb0 * rb1 * rb2); hnrb->Scale( scale_factor ); return hnrb; } -inline THn* rebin_visibility_map(const THn* hn, const int rbx, const int rby, const int rbz) { - const Int_t rbidx[3] = {rbx, rby, rbz}; +/** + * @brief Rebin the visibility map averaging the bin content after rebin + * + * @param hn visibility map + * @param rb0 bin group for axis 0 + * @param rb1 bin group for axis 1 + * @param rb2 bin group for axis 2 + * @return rebinned visibility map + */ +inline THn* rebin_visibility_map(const THn* hn, const int rb0, const int rb1, const int rb2) { + const Int_t rbidx[3] = {rb0, rb1, rb2}; THn* hnrb = hn->Rebin(rbidx); - float scale_factor = 1.0 / (rbx * rby * rbz); + float scale_factor = 1.0 / (rb0 * rb1 * rb2); hnrb->Scale( scale_factor ); return hnrb; } +/** + * @brief Get axis along which averaging the visibility (internal use) + */ inline const int get_projection_axis(const int& axis0, const int& axis1) { assert(axis0 != axis1); assert(axis1 < 2); @@ -55,7 +72,17 @@ inline const int get_projection_axis(const int& axis0, const int& axis1) { return (3-axis0-axis1); } -inline TH2D* draw_projection(const THnBase* hn, +/** + * @brief Compute a 2D average visibility profile + * + * @param hn visibility map + * @param axis0 axis index on the projection y-axis + * @param axis1 axis index on the projection x-axis + * @param x0 lower edge of the projected axis range + * @param x1 upper edge of the projected axis range + * @return projected visibility + */ +inline TH2D* get_projection(const THnBase* hn, const int& axis0, const int& axis1, const double& x0, @@ -74,76 +101,5 @@ inline TH2D* draw_projection(const THnBase* hn, return h_proj; } -inline void CanvasPartition(TCanvas *C,const Int_t Nx,const Int_t Ny, - Float_t lMargin, Float_t rMargin, - Float_t bMargin, Float_t tMargin) -{ - if (!C) return; - // Setup Pad layout: - Float_t vSpacing = 0.0; - Float_t vStep = (1.- bMargin - tMargin - (Ny-1) * vSpacing) / Ny; - Float_t hSpacing = 0.0; - Float_t hStep = (1.- lMargin - rMargin - (Nx-1) * hSpacing) / Nx; - Float_t vposd,vposu,vmard,vmaru,vfactor; - Float_t hposl,hposr,hmarl,hmarr,hfactor; - for (Int_t i=0;icd(0); - char name[16]; - sprintf(name,"pad_%i_%i",i,j); - TPad *pad = (TPad*) gROOT->FindObject(name); - if (pad) delete pad; - pad = new TPad(name,"",hposl,vposd,hposr,vposu); - pad->SetLeftMargin(hmarl); - pad->SetRightMargin(hmarr); - pad->SetBottomMargin(vmard); - pad->SetTopMargin(vmaru); - pad->SetFrameBorderMode(0); - pad->SetBorderMode(0); - pad->SetBorderSize(0); - pad->Draw(); - } - } -} - - #endif /* end of include guard DUNEVISUTILS_HPP */ diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C index 2874400d..e7b3dd80 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C @@ -2,6 +2,7 @@ * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) * @file : ttree_to_th3.C * @created : 2024-10-18 08:13 + * @description : Fill 3D visibility maps from TTree produced by the PhotonVisibilityExport */ #include @@ -31,7 +32,21 @@ #include "TRandom3.h" #include "TROOT.h" -void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_t vis_threshold = 1e-4) +/** + * Create THn visibility maps from the PhotonVisibilityExport_module's output + * + * Create visibility maps for the entire detector volume and for each individual OpDet. + * To reduce the resources needed for handling the visibility of the full set of + * OpDet, we use three-dimensional THnSparse to store the visiblity map, rounding + * the visibility to 0.0 when below the threshold set in vis_threshold. + * + * @param input_file input file path + * @param visexport_label label used for the PhotonVisibilityExport module + * @param vis_threshold Set visibility to 0 if below this value. (Only for single OpDet maps) + */ +void ttree_to_th3(const TString input_file, + const TString visexport_label, + const Double_t vis_threshold = 1e-4) { gStyle->SetTitleSize(0.06, "xyz"); gStyle->SetLabelSize(0.06, "xyz"); @@ -40,7 +55,7 @@ void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_ TTree* tOpDet = nullptr; TTree* tVisTPC = nullptr; TTree* tVisBuf = nullptr; - TDirectoryFile* dunevis_dir = file->Get(vis_dir); + TDirectoryFile* dunevis_dir = file->Get(visexport_label); dunevis_dir->cd(); TH1D* h_tmp[3] = {0}; @@ -65,9 +80,9 @@ void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_ TTreeReader reader(tVisTPC); TTreeReaderArray _xyz(reader, "coords"); TTreeReaderArray _opdet_visd(reader, "opDet_visDirect"); - //TTreeReaderArray _opdet_visr(reader, "opDet_visReflct"); + TTreeReaderArray _opdet_visr(reader, "opDet_visReflct"); TTreeReaderValue _vis_d(reader, "total_visDirect"); - //TTreeReaderValue _vis_r(reader, "total_visReflct"); + TTreeReaderValue _vis_r(reader, "total_visReflct"); const Long64_t N_OPDET = tOpDet->GetEntries(); @@ -99,9 +114,9 @@ void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_ TTreeReader readerBuff(tVisBuf); TTreeReaderArray _buf_xyz(readerBuff, "coords"); TTreeReaderArray _buf_opdet_visd(readerBuff, "opDet_visDirectBuff"); - //TTreeReaderArray _buf_opdet_visr(readerBuff, "opDet_visReflctBuff"); + TTreeReaderArray _buf_opdet_visr(readerBuff, "opDet_visReflctBuff"); TTreeReaderValue _buf_vis_d(readerBuff, "total_visDirectBuff"); - //TTreeReaderValue _buf_vis_r(readerBuff, "total_visReflctBuff"); + TTreeReaderValue _buf_vis_r(readerBuff, "total_visReflctBuff"); while( readerBuff.Next() ) { double xyz[3] = {_buf_xyz[0], _buf_xyz[1], _buf_xyz[2]}; @@ -118,7 +133,7 @@ void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_ TString output_name = input_file; output_name.Resize( output_name.Index(".root") ); - output_name += "_"+vis_dir+".root"; + output_name += "_"+visexport_label+".root"; TFile* output_h3 = new TFile(output_name, "recreate"); h3total_vis->Write(); @@ -130,15 +145,3 @@ void ttree_to_th3(const TString input_file, const TString vis_dir, const Double_ return; } -TH3F* rebin_visibility_map(const TH3F* h3, const int rbx, const int rby, const int rbz) -{ - TH3F* h3rb = (TH3F*)h3->Clone(); - h3rb = (TH3F*)h3rb->Rebin3D(rbx, rby, rbz, Form("%s_rb_%i_%i_%i", h3->GetName(), rbx, rby, rbz)); - - float scale_factor = 1.0 / (rbx * rby *rbz); - h3rb->Scale( scale_factor ); - - return h3rb; -} - - From d5b96c58ccafffb2487e80e31a841560eb4bae37 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Wed, 23 Oct 2024 07:29:13 -0500 Subject: [PATCH 4/9] better handling of reflected light --- duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl | 2 +- .../PhotonPropagation/PhotonVisibilityExport_module.cc | 7 ++++++- 2 files changed, 7 insertions(+), 2 deletions(-) diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl index 450c714d..3681b49f 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl @@ -13,7 +13,7 @@ standard_photovis: vis_model: "semianalytical" # visibility model to be used # (pick one between "semianalytical", "compgraph") do_refl: @local::dunefd_pdfastsim_par_ar_refl.DoReflectedLight - do_include_anode_refl: "false" + do_include_anode_refl: false vuvhitspars: @local::dunefd_pdfastsim_par_ar.VUVHits vishitspars: @local::dunefd_pdfastsim_par_ar_refl.VISHits tfloaderpars: @local::dunevd_pdfastsim_ann_ar.TFLoaderTool diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc index f8e1100a..f4d76bd3 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc @@ -159,6 +159,9 @@ namespace opdet { // create the photo-detector visibility model if (kVisModel == kSemiAnalytical) { printf("Creating Semi-analytical visibility model\n"); + (fDoReflectedLight) ? + printf("Reflections included\n") : + printf("Reflections NOT included\n"); fVisibilityModel = std::make_unique( fVUVHitsParams, fVISHitsParams, fDoReflectedLight, fIncludeAnodeReflections, fUseXeAbsorption @@ -436,7 +439,9 @@ namespace opdet { if (kVisModel == kSemiAnalytical ) { fVisibilityModel->detectedDirectVisibilities (opdetvis_dir, xspot); - fVisibilityModel->detectedReflectedVisibilities(opdetvis_rfl, xspot, fIncludeAnodeReflections); + if (fDoReflectedLight) { + fVisibilityModel->detectedReflectedVisibilities(opdetvis_rfl, xspot, fIncludeAnodeReflections); + } } else if (kVisModel == kCompGraph ) { std::vector pos_tmp {xspot.x(), xspot.y(), xspot.z()}; From cf61952f0b486bb053958e34559fca5583f9f093 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Wed, 23 Oct 2024 07:45:39 -0500 Subject: [PATCH 5/9] split filling of vis map to reduce memory impact --- .../VisibilityMapTools/ttree_to_th3.C | 178 +++++++++++------- 1 file changed, 111 insertions(+), 67 deletions(-) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C index e7b3dd80..c42deaf4 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C @@ -7,30 +7,22 @@ #include #include -#include #include -#include - -#include "RtypesCore.h" -#include "TFile.h" -#include "TTree.h" -#include "TTreeReader.h" -#include "TTreeReaderArray.h" -#include "TTreeReaderValue.h" -#include "TDirectoryFile.h" -#include "TSystem.h" -#include "TStyle.h" -#include "TCanvas.h" -#include "TClonesArray.h" -#include "TH1D.h" -#include "TH3F.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include #include #include -#include "TH2D.h" -#include "TGraphErrors.h" -#include "TF1.h" -#include "TRandom3.h" -#include "TROOT.h" +#include +#include /** * Create THn visibility maps from the PhotonVisibilityExport_module's output @@ -43,10 +35,12 @@ * @param input_file input file path * @param visexport_label label used for the PhotonVisibilityExport module * @param vis_threshold Set visibility to 0 if below this value. (Only for single OpDet maps) + * @param split_fill Split the filling of OpDet visivility in N runs to reduce memory impact (default 5) */ void ttree_to_th3(const TString input_file, const TString visexport_label, - const Double_t vis_threshold = 1e-4) + const Double_t vis_threshold = 1e-4, + const UInt_t split_fill = 5) { gStyle->SetTitleSize(0.06, "xyz"); gStyle->SetLabelSize(0.06, "xyz"); @@ -58,11 +52,11 @@ void ttree_to_th3(const TString input_file, TDirectoryFile* dunevis_dir = file->Get(visexport_label); dunevis_dir->cd(); + const TString axis_label[3] = {"x [cm]", "y [cm]", "z [cm]"}; TH1D* h_tmp[3] = {0}; h_tmp[0] = dunevis_dir->Get("hgrid0"); h_tmp[1] = dunevis_dir->Get("hgrid1"); h_tmp[2] = dunevis_dir->Get("hgrid2"); - const int nbins[3] = {h_tmp[0]->GetNbinsX(), h_tmp[1]->GetNbinsX(), h_tmp[2]->GetNbinsX()}; const double xmin[3] ={ h_tmp[0]->GetXaxis()->GetXmin(), @@ -84,63 +78,113 @@ void ttree_to_th3(const TString input_file, TTreeReaderValue _vis_d(reader, "total_visDirect"); TTreeReaderValue _vis_r(reader, "total_visReflct"); - const Long64_t N_OPDET = tOpDet->GetEntries(); - - THnF* h3total_vis = new THnF("h3VisMap", "Visibility Map", 3, nbins, xmin, xmax); - std::vector h3opDet; - h3opDet.reserve(N_OPDET); - for (size_t iopdet = 0; iopdet < N_OPDET; iopdet++) { - TString name = Form("h3VisMap_opDet%ld", iopdet); - TString titl = Form("Visibility Map - OpDet %ld", iopdet); - h3opDet.emplace_back(new THnSparseF(name, titl, 3, nbins, xmin, xmax)); + const size_t N_OPDET = tOpDet->GetEntries(); + size_t nOpDetPerRun = N_OPDET; + UInt_t nRuns = 1; + if (split_fill > 1) { + nRuns = split_fill; + nOpDetPerRun = std::ceil( N_OPDET / (double)nRuns); + } + std::vector> OpDetChunks(nRuns); + size_t iopdet = 0; + for (size_t irun = 0; irun < nRuns; irun++) { + std::vector& opdets_ = OpDetChunks.at(irun); + opdets_.reserve(nOpDetPerRun); + while (iopdet < N_OPDET && opdets_.size() <= nOpDetPerRun) { + opdets_.emplace_back(iopdet); + iopdet++; + } } - TH1D* h_opdet_vis_tpc = new TH1D("h_opdet_vis_tpc", "opdet visibility (TPC)", 100, -15, 1); - TH1D* h_opdet_vis_buf = new TH1D("h_opdet_vis_buf", "opdet visibility (buffer)", 100, -15, 1); + TString output_name = input_file; + output_name.Resize( output_name.Index(".root") ); + output_name += "_"+visexport_label+".root"; + TFile* output_h3 = new TFile(output_name, "recreate"); + + size_t irun = 0; + for (const auto& opdets : OpDetChunks) { + printf("Processing OpDet chunk %ld - [%ld - %ld]\n", irun, opdets.front(), opdets.back()); + + TClonesArray h3opDet("THnSparseF", opdets.size()); + size_t iopdet_idx = 0; + for (const auto& iopdet : opdets) { + TString name = Form("h3VisMap_opDet%ld", iopdet); + TString titl = Form("Visibility Map - OpDet %ld", iopdet); + new (h3opDet[iopdet_idx]) THnSparseF(name, titl, 3, nbins, xmin, xmax); + for (size_t idim = 0; idim < 3; idim++) { + ((THnSparseF*)h3opDet[iopdet_idx])->GetAxis(idim)->SetTitle(axis_label[idim]); + } + iopdet_idx++; + } - while( reader.Next() ) { - double xyz[3] = {_xyz[0], _xyz[1], _xyz[2]}; - h3total_vis->Fill(xyz, *(_vis_d)); - for (Long64_t i = 0; i < N_OPDET; i++) { - h_opdet_vis_tpc->Fill( log10(_opdet_visd[i]) ); - if (_opdet_visd[i] > vis_threshold) { - h3opDet[i]->Fill(xyz, _opdet_visd[i]); + THnF* hnTotal = nullptr; + if (irun == 0) { + hnTotal = new THnF("h3VisMap", "Detector visibility map", 3, nbins, xmin, xmax); + for (size_t idim = 0; idim < 3; idim++) { + hnTotal->GetAxis(idim)->SetTitle(axis_label[idim]); } } - } - if (tVisBuf) { - if (tVisBuf->GetEntries() > 0) { - TTreeReader readerBuff(tVisBuf); - TTreeReaderArray _buf_xyz(readerBuff, "coords"); - TTreeReaderArray _buf_opdet_visd(readerBuff, "opDet_visDirectBuff"); - TTreeReaderArray _buf_opdet_visr(readerBuff, "opDet_visReflctBuff"); - TTreeReaderValue _buf_vis_d(readerBuff, "total_visDirectBuff"); - TTreeReaderValue _buf_vis_r(readerBuff, "total_visReflctBuff"); - - while( readerBuff.Next() ) { - double xyz[3] = {_buf_xyz[0], _buf_xyz[1], _buf_xyz[2]}; - h3total_vis->Fill(xyz, *(_buf_vis_d)); - for (Long64_t i = 0; i < N_OPDET; i++) { - h_opdet_vis_buf->Fill( log10(_buf_opdet_visd[i]) ); - if ( _buf_opdet_visd[i] > vis_threshold ) { - h3opDet[i]->Fill(xyz, _buf_opdet_visd[i]); + reader.Restart(); + + while (reader.Next()) { + double xyz[3] = {_xyz[0], _xyz[1], _xyz[2]}; + iopdet_idx = 0; + for (const auto& iopdet : opdets) { + double& vis = _opdet_visd[iopdet]; + if (vis > vis_threshold) { + ((THnSparseF*)h3opDet[iopdet_idx])->Fill(xyz, vis); + } + iopdet_idx++; + } + + if (irun == 0) { + hnTotal->Fill( xyz, *(_vis_d) ); + } + } + + if (tVisBuf) { + if (tVisBuf->GetEntries() > 0) { + TTreeReader readerBuff(tVisBuf); + TTreeReaderArray _buf_xyz(readerBuff, "coords"); + TTreeReaderArray _buf_opdet_visd(readerBuff, "opDet_visDirectBuff"); + TTreeReaderArray _buf_opdet_visr(readerBuff, "opDet_visReflctBuff"); + TTreeReaderValue _buf_vis_d(readerBuff, "total_visDirectBuff"); + TTreeReaderValue _buf_vis_r(readerBuff, "total_visReflctBuff"); + + while( readerBuff.Next() ) { + double xyz[3] = {_buf_xyz[0], _buf_xyz[1], _buf_xyz[2]}; + iopdet_idx = 0; + for (const auto& iopdet : opdets) { + double& vis = _opdet_visd[iopdet]; + if ( vis > vis_threshold ) { + ((THnSparseF*)h3opDet[iopdet_idx])->Fill(xyz, vis); + } + iopdet_idx++; } + + if (irun == 0) { + hnTotal->Fill( xyz, *(_buf_vis_d) ); + } + } } } - } - TString output_name = input_file; - output_name.Resize( output_name.Index(".root") ); - output_name += "_"+visexport_label+".root"; + output_h3->cd(); + if (irun == 0) { + hnTotal->Write(); + delete hnTotal; + } + for (const auto& h3 : h3opDet) { + h3->Write(); + } + h3opDet.Clear(); - TFile* output_h3 = new TFile(output_name, "recreate"); - h3total_vis->Write(); - for (auto& h3 : h3opDet) { - h3->Write(); + irun++; } - output_h3->Close(); + + file->Close(); return; } From 4a76c13045f173f724f0c727f8cf8eff981d6291 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Wed, 23 Oct 2024 07:47:50 -0500 Subject: [PATCH 6/9] Add README and example fhicls --- .../VisibilityMapTools/README.md | 20 ++++++++ .../fcl/lightmap_dune_fdhd_1x2x6.fcl | 51 +++++++++++++++++++ .../fcl/lightmap_dune_fdvd_1x8x14_ar.fcl | 51 +++++++++++++++++++ .../fcl/lightmap_dune_fdvd_1x8x14_xe.fcl | 51 +++++++++++++++++++ 4 files changed, 173 insertions(+) create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/README.md create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl create mode 100644 duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/README.md b/duneopdet/PhotonPropagation/VisibilityMapTools/README.md new file mode 100644 index 00000000..1015db4b --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/README.md @@ -0,0 +1,20 @@ +--- +author: Daniele Guffanti - University & INFN Milano-Bicocca (daniele.guffanti@mib.infn.it) +date: 2024-10-22 08:25 +--- +# VisibilityMap Tools + +This folder contains a few example macros for creating 3D/2D visibility maps +using the `PhotonVisibilityExport_module`. + +Running the macro `ttree_to_th3.C` will produce a file containing the detector visibility +map stored in a three-dimensional THnF and the individual visibility maps for each OpDet +as `THnSparseF` objects. + +The file `DUNEVisUtils.hpp` containes a few useful functions to handle these histograms, +in particular to "rebin" the maps for a better display and to profile them along a given +dimension. + +Example fhicl configurations can be found in the `fcl` folder. Note that in case of +Xe doping, the visivility must be exported for Ar and Xe separately and combined +in a later stage with the appropriate weights. diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl new file mode 100644 index 00000000..8d6b6b0d --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl @@ -0,0 +1,51 @@ +#include "PDFastSim_dune.fcl" +#include "services_dunefd_horizdrift_1x2x6.fcl" +#include "PhotonVisibilityExport.fcl" + +process_name: Ana + +services: +{ + # Load the service that manages root files for histograms. + TFileService: { fileName: "dunevis_fdhd_1x2x6_test.root" } + TimeTracker: {} + RandomNumberGenerator: {} + MemoryTracker: { } # default is one + message: @local::dune_message_services_prod + + ### @table::dunefd_refactored_services + @table::dunefd_1x2x6_services + PhotonVisibilityService: @local::dune10kt_1x2x6_photonvisibilityservice +} + +### Use the 1x2x6 geometry ### +services.Geometry: @local::dune10kt_1x2x6_refactored_geo + +#source is now a root file +source: +{ + module_type: RootInput + fileNames : ["dunevis_hd_1x2x6_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisAr: @local::standard_photovis + } + + anapath: [ photovisAr ] + end_paths: [ anapath ] +} + +physics.analyzers.photovisAr.vis_model: "semianalytical" +physics.analyzers.photovisAr.vuvhitspars: @local::dunefd_pdfastsim.VUVHits +physics.analyzers.photovisAr.do_include_buffer: "true" +physics.analyzers.photovisAr.voxel_dx: "15.0" # voxel x mesh step (in cm) +physics.analyzers.photovisAr.voxel_dy: "15.0" # voxel y mesh step (in cm) +physics.analyzers.photovisAr.voxel_dz: "15.0" # voxel z mesh step (in cm) + + diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl new file mode 100644 index 00000000..4c64d4dc --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl @@ -0,0 +1,51 @@ +#include "PDFastSim_dune.fcl" +#include "services_dunefd_vertdrift_1x8x14.fcl" +#include "PhotonVisibilityExport.fcl" + +process_name: Ana + +services: +{ + # Load the service that manages root files for histograms. + TFileService: { fileName: "dunevis_fdvd_1x8x14_ar.root" } + TimeTracker: {} + RandomNumberGenerator: {} + MemoryTracker: { } # default is one + message: @local::dune_message_services_prod + + ### @table::dunefd_refactored_services + @table::dunefdvd_1x8x14_3view_30deg_simulation_services +} + +### Use the 1x2x6 geometry ### +services.Geometry: @local::dunevd10kt_1x8x14_3view_30deg_v5_geo + +#source is now a root file +source: +{ + module_type: RootInput + fileNames : ["dunevis_vd_1x8x14_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisAr: @local::standard_photovis + } + + anapath: [ photovisAr ] + end_paths: [ anapath ] +} + +physics.analyzers.photovisAr.vis_model: "semianalytical" +physics.analyzers.photovisAr.vuvhitspars: @local::dunevd_pdfastsim_par_ar.VUVHits +physics.analyzers.photovisAr.do_refl: @local::dunevd_pdfastsim_par_ar.DoReflectedLight +physics.analyzers.photovisAr.do_include_buffer: true +physics.analyzers.photovisAr.voxel_dx: 15.0 # voxel x mesh step (in cm) +physics.analyzers.photovisAr.voxel_dy: 15.0 # voxel y mesh step (in cm) +physics.analyzers.photovisAr.voxel_dz: 15.0 # voxel z mesh step (in cm) + + diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl new file mode 100644 index 00000000..7f1e01fe --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl @@ -0,0 +1,51 @@ +#include "PDFastSim_dune.fcl" +#include "services_dunefd_vertdrift_1x8x14.fcl" +#include "PhotonVisibilityExport.fcl" + +process_name: Ana + +services: +{ + # Load the service that manages root files for histograms. + TFileService: { fileName: "dunevis_fdvd_1x8x14_xe.root" } + TimeTracker: {} + RandomNumberGenerator: {} + MemoryTracker: { } # default is one + message: @local::dune_message_services_prod + + ### @table::dunefd_refactored_services + @table::dunefdvd_1x8x14_3view_30deg_simulation_services +} + +### Use the 1x2x6 geometry ### +services.Geometry: @local::dunevd10kt_1x8x14_3view_30deg_v5_geo +services.PhotonVisibilityService: @local::dune10kt_vd_photonvisibilityservice_Xe # changing to the xenon libaray. + +#source is now a root file +source: +{ + module_type: RootInput + fileNames : ["dunevis_vd_1x8x14_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisXe: @local::standard_photovis + } + + anapath: [ photovisXe ] + end_paths: [ anapath ] +} + +physics.analyzers.photovisXe.vis_model: "semianalytical" +physics.analyzers.photovisXe.do_include_buffer: "true" +physics.analyzers.photovisXe.vuvhitspars: @local::dunevd_pdfastsim_par_xe.VUVHits +physics.analyzers.photovisXe.do_refl: false +physics.analyzers.photovisXe.do_include_anode_refl: false +physics.analyzers.photovisXe.voxel_dx: "15.0" # voxel x mesh step (in cm) +physics.analyzers.photovisXe.voxel_dy: "15.0" # voxel y mesh step (in cm) +physics.analyzers.photovisXe.voxel_dz: "15.0" # voxel z mesh step (in cm) From 2a25469caf0233ec9a72266117ddf85d93803e67 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Fri, 8 Nov 2024 04:58:12 -0600 Subject: [PATCH 7/9] Add FD HD/VD configuration in standard fhicl --- .../PhotonVisibilityExport.fcl | 32 +++++++++++++------ .../fcl/lightmap_dune_fdhd_1x2x6.fcl | 11 +++---- .../fcl/lightmap_dune_fdvd_1x8x14_ar.fcl | 15 ++++----- .../fcl/lightmap_dune_fdvd_1x8x14_xe.fcl | 18 +++++------ 4 files changed, 43 insertions(+), 33 deletions(-) diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl index 3681b49f..ba3a6c0d 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl @@ -2,24 +2,36 @@ BEGIN_PROLOG -standard_photovis: +standard_photovisexport: { - module_type: "PhotonVisibilityExport" - photovis_label: "photovis" - photovis_map: "photovis" - voxel_dx: "10.00" # voxel x mesh step (in cm) - voxel_dy: "10.00" # voxel y mesh step (in cm) - voxel_dz: "10.00" # voxel z mesh step (in cm) - vis_model: "semianalytical" # visibility model to be used - # (pick one between "semianalytical", "compgraph") + module_type: "PhotonVisibilityExport" + photovis_label: "photovisexport" + photovis_map: "photovisexport" + voxel_dx: "10.00" # voxel x mesh step (in cm) + voxel_dy: "10.00" # voxel y mesh step (in cm) + voxel_dz: "10.00" # voxel z mesh step (in cm) + vis_model: "semianalytical" # visibility model to be used + # (pick one between "semianalytical", "compgraph") do_refl: @local::dunefd_pdfastsim_par_ar_refl.DoReflectedLight do_include_anode_refl: false vuvhitspars: @local::dunefd_pdfastsim_par_ar.VUVHits vishitspars: @local::dunefd_pdfastsim_par_ar_refl.VISHits tfloaderpars: @local::dunevd_pdfastsim_ann_ar.TFLoaderTool - do_include_buffer: false + do_include_buffer: true # include buffer region using the PhotonVisibilityService } +photovisexport_fdhd: @local::standard_photovisexport + +photovisexport_fdvd_ar: @local::standard_photovisexport +photovisexport_fdvd_ar.vuvhitspars: @local::dunevd_pdfastsim_par_ar.VUVHits +photovisexport_fdvd_ar.do_refl: @local::dunevd_pdfastsim_par_ar.DoReflectedLight +photovisexport_fdvd_ar.do_include_anode_refl: @local::dunevd_pdfastsim_par_ar.IncludeAnodeReflections + +photovisexport_fdvd_xe: @local::standard_photovisexport +photovisexport_fdvd_xe.vuvhitspars: @local::dunevd_pdfastsim_par_xe.VUVHits +photovisexport_fdvd_xe.do_refl: @local::dunevd_pdfastsim_par_xe.DoReflectedLight +photovisexport_fdvd_xe.do_include_anode_refl: @local::dunevd_pdfastsim_par_xe.IncludeAnodeReflections + END_PROLOG diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl index 8d6b6b0d..59fcf86b 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl @@ -13,7 +13,6 @@ services: MemoryTracker: { } # default is one message: @local::dune_message_services_prod - ### @table::dunefd_refactored_services @table::dunefd_1x2x6_services PhotonVisibilityService: @local::dune10kt_1x2x6_photonvisibilityservice } @@ -21,7 +20,6 @@ services: ### Use the 1x2x6 geometry ### services.Geometry: @local::dune10kt_1x2x6_refactored_geo -#source is now a root file source: { module_type: RootInput @@ -38,12 +36,13 @@ physics: } anapath: [ photovisAr ] - end_paths: [ anapath ] + end_paths: [ anapath ] } -physics.analyzers.photovisAr.vis_model: "semianalytical" -physics.analyzers.photovisAr.vuvhitspars: @local::dunefd_pdfastsim.VUVHits -physics.analyzers.photovisAr.do_include_buffer: "true" +# physics.analyzers.photovisAr.vis_model: "semianalytical" +# physics.analyzers.photovisAr.vuvhitspars: @local::dunefd_pdfastsim.VUVHits +# physics.analyzers.photovisAr.do_include_buffer: "true" +physics.analyzers.photovisAr: @local::photovisexport_fdhd physics.analyzers.photovisAr.voxel_dx: "15.0" # voxel x mesh step (in cm) physics.analyzers.photovisAr.voxel_dy: "15.0" # voxel y mesh step (in cm) physics.analyzers.photovisAr.voxel_dz: "15.0" # voxel z mesh step (in cm) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl index 4c64d4dc..971fdcb9 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl @@ -13,14 +13,12 @@ services: MemoryTracker: { } # default is one message: @local::dune_message_services_prod - ### @table::dunefd_refactored_services @table::dunefdvd_1x8x14_3view_30deg_simulation_services } -### Use the 1x2x6 geometry ### +### Use the appropriate geometry ### services.Geometry: @local::dunevd10kt_1x8x14_3view_30deg_v5_geo -#source is now a root file source: { module_type: RootInput @@ -37,13 +35,14 @@ physics: } anapath: [ photovisAr ] - end_paths: [ anapath ] + end_paths: [ anapath ] } -physics.analyzers.photovisAr.vis_model: "semianalytical" -physics.analyzers.photovisAr.vuvhitspars: @local::dunevd_pdfastsim_par_ar.VUVHits -physics.analyzers.photovisAr.do_refl: @local::dunevd_pdfastsim_par_ar.DoReflectedLight -physics.analyzers.photovisAr.do_include_buffer: true +#physics.analyzers.photovisAr.vis_model: "semianalytical" +#physics.analyzers.photovisAr.vuvhitspars: @local::dunevd_pdfastsim_par_ar.VUVHits +#physics.analyzers.photovisAr.do_refl: @local::dunevd_pdfastsim_par_ar.DoReflectedLight +#physics.analyzers.photovisAr.do_include_buffer: true +physics.analyzers.photovisAr: @local::photovisexport_fdvd_ar physics.analyzers.photovisAr.voxel_dx: 15.0 # voxel x mesh step (in cm) physics.analyzers.photovisAr.voxel_dy: 15.0 # voxel y mesh step (in cm) physics.analyzers.photovisAr.voxel_dz: 15.0 # voxel z mesh step (in cm) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl index 7f1e01fe..33080b56 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl @@ -13,15 +13,14 @@ services: MemoryTracker: { } # default is one message: @local::dune_message_services_prod - ### @table::dunefd_refactored_services @table::dunefdvd_1x8x14_3view_30deg_simulation_services } -### Use the 1x2x6 geometry ### +### Use the appropriate geometry ### services.Geometry: @local::dunevd10kt_1x8x14_3view_30deg_v5_geo +### Setup the visibility service for Xe services.PhotonVisibilityService: @local::dune10kt_vd_photonvisibilityservice_Xe # changing to the xenon libaray. -#source is now a root file source: { module_type: RootInput @@ -38,14 +37,15 @@ physics: } anapath: [ photovisXe ] - end_paths: [ anapath ] + end_paths: [ anapath ] } -physics.analyzers.photovisXe.vis_model: "semianalytical" -physics.analyzers.photovisXe.do_include_buffer: "true" -physics.analyzers.photovisXe.vuvhitspars: @local::dunevd_pdfastsim_par_xe.VUVHits -physics.analyzers.photovisXe.do_refl: false -physics.analyzers.photovisXe.do_include_anode_refl: false +# physics.analyzers.photovisXe.vis_model: "semianalytical" +# physics.analyzers.photovisXe.do_include_buffer: "true" +# physics.analyzers.photovisXe.vuvhitspars: @local::dunevd_pdfastsim_par_xe.VUVHits +# physics.analyzers.photovisXe.do_refl: false +# physics.analyzers.photovisXe.do_include_anode_refl: false +physics.analyzers.photovisXe: @local::photovisexport_fdvd_xe physics.analyzers.photovisXe.voxel_dx: "15.0" # voxel x mesh step (in cm) physics.analyzers.photovisXe.voxel_dy: "15.0" # voxel y mesh step (in cm) physics.analyzers.photovisXe.voxel_dz: "15.0" # voxel z mesh step (in cm) From f0dc76a1c6629bcf24b3982f73392f4b1c2a8b17 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Fri, 8 Nov 2024 10:45:53 -0600 Subject: [PATCH 8/9] fixed typo in standard configuration sourcing --- .../VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl | 6 +----- .../fcl/lightmap_dune_fdvd_1x8x14_ar.fcl | 7 +------ .../fcl/lightmap_dune_fdvd_1x8x14_xe.fcl | 8 +------- 3 files changed, 3 insertions(+), 18 deletions(-) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl index 59fcf86b..3132d08e 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl @@ -32,17 +32,13 @@ physics: { analyzers: { - photovisAr: @local::standard_photovis + photovisAr: @local::photovisexport_fdhd } anapath: [ photovisAr ] end_paths: [ anapath ] } -# physics.analyzers.photovisAr.vis_model: "semianalytical" -# physics.analyzers.photovisAr.vuvhitspars: @local::dunefd_pdfastsim.VUVHits -# physics.analyzers.photovisAr.do_include_buffer: "true" -physics.analyzers.photovisAr: @local::photovisexport_fdhd physics.analyzers.photovisAr.voxel_dx: "15.0" # voxel x mesh step (in cm) physics.analyzers.photovisAr.voxel_dy: "15.0" # voxel y mesh step (in cm) physics.analyzers.photovisAr.voxel_dz: "15.0" # voxel z mesh step (in cm) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl index 971fdcb9..2b3def50 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl @@ -31,18 +31,13 @@ physics: { analyzers: { - photovisAr: @local::standard_photovis + photovisAr: @local::photovisexport_fdvd_ar } anapath: [ photovisAr ] end_paths: [ anapath ] } -#physics.analyzers.photovisAr.vis_model: "semianalytical" -#physics.analyzers.photovisAr.vuvhitspars: @local::dunevd_pdfastsim_par_ar.VUVHits -#physics.analyzers.photovisAr.do_refl: @local::dunevd_pdfastsim_par_ar.DoReflectedLight -#physics.analyzers.photovisAr.do_include_buffer: true -physics.analyzers.photovisAr: @local::photovisexport_fdvd_ar physics.analyzers.photovisAr.voxel_dx: 15.0 # voxel x mesh step (in cm) physics.analyzers.photovisAr.voxel_dy: 15.0 # voxel y mesh step (in cm) physics.analyzers.photovisAr.voxel_dz: 15.0 # voxel z mesh step (in cm) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl index 33080b56..a8bfde3c 100644 --- a/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl @@ -33,19 +33,13 @@ physics: { analyzers: { - photovisXe: @local::standard_photovis + photovisXe: @local::photovisexport_fdvd_xe } anapath: [ photovisXe ] end_paths: [ anapath ] } -# physics.analyzers.photovisXe.vis_model: "semianalytical" -# physics.analyzers.photovisXe.do_include_buffer: "true" -# physics.analyzers.photovisXe.vuvhitspars: @local::dunevd_pdfastsim_par_xe.VUVHits -# physics.analyzers.photovisXe.do_refl: false -# physics.analyzers.photovisXe.do_include_anode_refl: false -physics.analyzers.photovisXe: @local::photovisexport_fdvd_xe physics.analyzers.photovisXe.voxel_dx: "15.0" # voxel x mesh step (in cm) physics.analyzers.photovisXe.voxel_dy: "15.0" # voxel y mesh step (in cm) physics.analyzers.photovisXe.voxel_dz: "15.0" # voxel z mesh step (in cm) From 2b3f2ca40ad7a8f1842fd3020142b54123f79e78 Mon Sep 17 00:00:00 2001 From: Daniele Guffanti Date: Thu, 14 Nov 2024 11:54:16 -0600 Subject: [PATCH 9/9] Add variable to set visibility sampling size --- .../PhotonVisibilityExport.fcl | 1 + .../PhotonVisibilityExport_module.cc | 82 ++++++++++++++----- 2 files changed, 62 insertions(+), 21 deletions(-) diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl index ba3a6c0d..201038cf 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl @@ -12,6 +12,7 @@ standard_photovisexport: voxel_dz: "10.00" # voxel z mesh step (in cm) vis_model: "semianalytical" # visibility model to be used # (pick one between "semianalytical", "compgraph") + n_vis_samplings: "5" # nr of visibility samplings inside the voxel (min 1) do_refl: @local::dunefd_pdfastsim_par_ar_refl.DoReflectedLight do_include_anode_refl: false vuvhitspars: @local::dunefd_pdfastsim_par_ar.VUVHits diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc index e29f92fd..dbd7533b 100644 --- a/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc @@ -96,6 +96,8 @@ namespace opdet { double fVoxelSizeY = {}; double fVoxelSizeZ = {}; + int fNSamplings = 1; + fhicl::ParameterSet fVUVHitsParams; fhicl::ParameterSet fVISHitsParams; fhicl::ParameterSet fTFLoaderPars; @@ -114,6 +116,10 @@ namespace opdet { void ExportVoxelGrid(); void ExportVisibility(); + const bool is_valid(const phot::MappedCounts_t& visibilities) const; + const bool is_valid(const std::vector& visibilities) const; + const geo::Point_t SampleVoxel(const geo::Point_t& vxc) const; + }; } // close opdet namespace @@ -146,13 +152,14 @@ namespace opdet { fTFLoaderPars = pset.get("tfloaderpars"); TString vis_model_str = pset.get("vis_model"); - if (vis_model_str.Contains("compgraph")) { kVisModel = kCompGraph; } else if (vis_model_str.Contains("semianalytical")) { kVisModel = kSemiAnalytical; } + + fNSamplings = pset.get("n_vis_samplings"); } void PhotonVisibilityExport::beginJob() { @@ -326,6 +333,32 @@ namespace opdet { return; } + const geo::Point_t PhotonVisibilityExport::SampleVoxel(const geo::Point_t& vxc) const + { + const geo::Vector_t delta( + gRandom->Uniform(-0.5*fVoxelSizeX, 0.5*fVoxelSizeX), + gRandom->Uniform(-0.5*fVoxelSizeY, 0.5*fVoxelSizeY), + gRandom->Uniform(-0.5*fVoxelSizeZ, 0.5*fVoxelSizeZ) ); + + const geo::Point_t xspot = vxc + delta; + + return xspot; + } + + const bool PhotonVisibilityExport::is_valid(const phot::MappedCounts_t& visibilities) const { + for (const auto& v : visibilities) { + if (v > 1.0) return false; + } + return true; + } + + const bool PhotonVisibilityExport::is_valid(const std::vector& visibilities) const { + for (const auto& v : visibilities) { + if (v > 1.0) return false; + } + return true; + } + void PhotonVisibilityExport::ExportVisibility() { const auto photonVisService = art::ServiceHandle(); @@ -364,7 +397,6 @@ namespace opdet { std::vector opdetvis_rfl; bool tpc_range_x = false, tpc_range_y = false, tpc_range_z = false; - const double voxelDim[3] = {fVoxelSizeX, fVoxelSizeY, fVoxelSizeZ}; // here we can loop over the points of the pre-defined grid double x_= 0, y_= 0, z_= 0; @@ -399,12 +431,27 @@ namespace opdet { } const geo::Point_t point( x_, y_, z_); - auto mapped_vis = photonVisService->GetAllVisibilities(point); - size_t iopdet = 0; - for (auto &vis : mapped_vis) { - opDet_visDirectBuff[iopdet] = vis; - total_visDirectBuff += vis; - iopdet++; + + int nvalid = 0; + while (nvalid < fNSamplings) { + const auto xpoint = SampleVoxel( point ); + auto mapped_vis = photonVisService->GetAllVisibilities(xpoint); + if ( is_valid(mapped_vis) == false ) { + continue; + } + + size_t iopdet = 0; + for (auto &vis : mapped_vis) { + opDet_visDirectBuff[iopdet] += vis; + total_visDirectBuff += vis; + iopdet++; + } + nvalid++; + } + + total_visDirectBuff = total_visDirectBuff / static_cast(fNSamplings); + for (size_t i = 0; i < fNOpDets; i++) { + opDet_visDirectBuff[i] /= static_cast(fNSamplings); } point.GetCoordinates( pointBuff_ ); @@ -420,15 +467,8 @@ namespace opdet { opDet_visReflct[i] = 0.0; } - const double n_samplings = 5.0; - - for (int i = 0; i < n_samplings; i++) { - const geo::Vector_t delta( - gRandom->Uniform(-0.5*voxelDim[0], 0.5*voxelDim[0]), - gRandom->Uniform(-0.5*voxelDim[1], 0.5*voxelDim[1]), - gRandom->Uniform(-0.5*voxelDim[2], 0.5*voxelDim[2]) ); - - const geo::Point_t xspot = vpoint + delta; + for (int i = 0; i < fNSamplings; i++) { + const geo::Point_t xspot = SampleVoxel( vpoint ); if (kVisModel == kSemiAnalytical ) { fVisibilityModel->detectedDirectVisibilities (opdetvis_dir, xspot); @@ -457,12 +497,12 @@ namespace opdet { } } - total_visDirect = total_visDirect / n_samplings; - total_visReflct = total_visReflct / n_samplings; + total_visDirect = total_visDirect / static_cast(fNSamplings); + total_visReflct = total_visReflct / static_cast(fNSamplings); for (size_t i = 0; i < fNOpDets; i++) { - opDet_visDirect[i] /= n_samplings; - opDet_visReflct[i] /= n_samplings; + opDet_visDirect[i] /= static_cast(fNSamplings); + opDet_visReflct[i] /= static_cast(fNSamplings); } // Fill tree