diff --git a/CMakeLists.txt b/CMakeLists.txt index 933db33f..7744367b 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -31,6 +31,9 @@ find_package( Boost REQUIRED ) find_package(ROOT REQUIRED) find_package( larevt REQUIRED EXPORT ) find_package( larsim REQUIRED EXPORT ) +find_package( eigen ) +find_package( tensorflow ) +find_package( larsimdnn ) find_package( larcore REQUIRED EXPORT ) find_package( lardata REQUIRED EXPORT ) find_package( lardataalg REQUIRED EXPORT ) @@ -42,6 +45,7 @@ find_package( dunecore REQUIRED EXPORT ) find_package( CLHEP REQUIRED ) find_package( Geant4 REQUIRED ) find_package( HDF5 REQUIRED ) + # macros for artdaq_dictionary and simple_plugin include(ArtDictionary) include(ArtMake) diff --git a/duneopdet/PhotonPropagation/CMakeLists.txt b/duneopdet/PhotonPropagation/CMakeLists.txt index f33fa51f..17e965e6 100644 --- a/duneopdet/PhotonPropagation/CMakeLists.txt +++ b/duneopdet/PhotonPropagation/CMakeLists.txt @@ -13,6 +13,7 @@ art_make( art::Framework_Services_Registry art_root_io::tfile_support ROOT::Core + ROOT::Tree art_root_io::TFileService_service art::Framework_Services_Optional_RandomNumberGenerator_service art::Persistency_Common @@ -28,8 +29,12 @@ art_make( larsim::LegacyLArG4 lardataobj::Simulation lardata::headers + larsim::PhotonPropagation_PhotonVisibilityService_service duneopdet::PhotonPropagation_PhotonVisibilityServiceS2_service larsim::Simulation + TensorFlow::cc + TensorFlow::framework + larsimdnn::TFLoaderTool nug4::ParticleNavigation larcorealg::Geometry larcore::Geometry_Geometry_service @@ -38,6 +43,7 @@ art_make( art::Framework_Services_Registry art_root_io::tfile_support ROOT::Core + ROOT::Tree art_root_io::TFileService_service art::Framework_Services_Optional_RandomNumberGenerator_service nurandom::RandomUtils_NuRandomService_service diff --git a/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl new file mode 100644 index 00000000..201038cf --- /dev/null +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport.fcl @@ -0,0 +1,38 @@ +#include "PDFastSim_dune.fcl" + +BEGIN_PROLOG + +standard_photovisexport: +{ + 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") + 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 + vishitspars: @local::dunefd_pdfastsim_par_ar_refl.VISHits + tfloaderpars: @local::dunevd_pdfastsim_ann_ar.TFLoaderTool + + 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/PhotonVisibilityExport_module.cc b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc new file mode 100644 index 00000000..dbd7533b --- /dev/null +++ b/duneopdet/PhotonPropagation/PhotonVisibilityExport_module.cc @@ -0,0 +1,525 @@ +/** + * @author : Daniele Guffanti (daniele.guffanti@mib.infn.it) + * @file : PhotonVisibilityExport_module.cc + * @created : 2024-10-22 07:50 + */ + +#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 + * @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: + 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 = {}; + + int fNSamplings = 1; + + 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(); + + 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 + +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; + } + + fNSamplings = pset.get("n_vis_samplings"); + } + + void PhotonVisibilityExport::beginJob() { + // 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 + ); + } + else if (kVisModel == kCompGraph) { + fTFGenerator = art::make_tool(fTFLoaderPars); + fTFGenerator->Initialization(); + } + + return; + } + + void PhotonVisibilityExport::analyze(const art::Event&) { + + if (fIsDone == false) { + 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]; + + 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 + UInt_t icryo = 0; + for (geo::CryostatGeo const& cryo : geom->Iterate()) { + printf("cryostat %u: [%g, %g, %g] - [%g, %g, %g]\n", icryo, + 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(); + icryo++; + } + + + // loop over all TPCs to get the active volume dimensions + for (geo::TPCGeo const& tpc : geom->Iterate()) { + auto point_center = tpc.GetCenter(); + 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]; + + 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"); + + // store info from Geometry service + fNOpDets = geom->NOpDets(); + for (size_t i : util::counter(fNOpDets)) { + geo::OpDetGeo const& opDet = geom->OpDetGeoFromOpDet(i); + auto center = opDet.GetCenter(); + center.GetCoordinates( opDetPos ); + opDetH = opDet.Height(); + opDetW = opDet.Width(); + opDetL = opDet.Length(); + + tOpDet->Fill(); + } + 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(); + // 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; + + // 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_); + + 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_ ); + 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; + } + + for (int i = 0; i < fNSamplings; i++) { + const geo::Point_t xspot = SampleVoxel( vpoint ); + + if (kVisModel == kSemiAnalytical ) { + fVisibilityModel->detectedDirectVisibilities (opdetvis_dir, xspot); + if (fDoReflectedLight) { + 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) { + 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 / static_cast(fNSamplings); + total_visReflct = total_visReflct / static_cast(fNSamplings); + + for (size_t i = 0; i < fNOpDets; i++) { + opDet_visDirect[i] /= static_cast(fNSamplings); + opDet_visReflct[i] /= static_cast(fNSamplings); + } + + // 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/DUNEVisUtils.hpp b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp new file mode 100644 index 00000000..a93bd52f --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/DUNEVisUtils.hpp @@ -0,0 +1,105 @@ +/** + * @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 + + +/** + * @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] = {rb0, rb1, rb2}; + THnSparse* hnrb = hn->Rebin(rbidx); + + float scale_factor = 1.0 / (rb0 * rb1 * rb2); + hnrb->Scale( scale_factor ); + return hnrb; +} + +/** + * @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 / (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); + assert(axis0 < 2); + + return (3-axis0-axis1); +} + +/** + * @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, + 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; +} + +#endif /* end of include guard DUNEVISUTILS_HPP */ + 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..3132d08e --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdhd_1x2x6.fcl @@ -0,0 +1,46 @@ +#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_1x2x6_services + PhotonVisibilityService: @local::dune10kt_1x2x6_photonvisibilityservice +} + +### Use the 1x2x6 geometry ### +services.Geometry: @local::dune10kt_1x2x6_refactored_geo + +source: +{ + module_type: RootInput + fileNames : ["dunevis_hd_1x2x6_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisAr: @local::photovisexport_fdhd + } + + anapath: [ photovisAr ] + end_paths: [ anapath ] +} + +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..2b3def50 --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_ar.fcl @@ -0,0 +1,45 @@ +#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::dunefdvd_1x8x14_3view_30deg_simulation_services +} + +### Use the appropriate geometry ### +services.Geometry: @local::dunevd10kt_1x8x14_3view_30deg_v5_geo + +source: +{ + module_type: RootInput + fileNames : ["dunevis_vd_1x8x14_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisAr: @local::photovisexport_fdvd_ar + } + + anapath: [ photovisAr ] + end_paths: [ anapath ] +} + +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..a8bfde3c --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/fcl/lightmap_dune_fdvd_1x8x14_xe.fcl @@ -0,0 +1,45 @@ +#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::dunefdvd_1x8x14_3view_30deg_simulation_services +} + +### 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: +{ + module_type: RootInput + fileNames : ["dunevis_vd_1x8x14_source.root"] + maxEvents: 1 +} + + +physics: +{ + analyzers: + { + photovisXe: @local::photovisexport_fdvd_xe + } + + anapath: [ photovisXe ] + end_paths: [ anapath ] +} + +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) diff --git a/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C new file mode 100644 index 00000000..c42deaf4 --- /dev/null +++ b/duneopdet/PhotonPropagation/VisibilityMapTools/ttree_to_th3.C @@ -0,0 +1,191 @@ +/** + * @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 +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +/** + * 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) + * @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 UInt_t split_fill = 5) +{ + 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(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(), + 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 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++; + } + } + + 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++; + } + + 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]); + } + } + + 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) ); + } + + } + } + } + + output_h3->cd(); + if (irun == 0) { + hnTotal->Write(); + delete hnTotal; + } + for (const auto& h3 : h3opDet) { + h3->Write(); + } + h3opDet.Clear(); + + irun++; + } + + file->Close(); + + return; +} +