diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ec0b2971..5298eda7 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -172,10 +172,10 @@ jobs: echo "[CI] RUN SIMULATION:" simulate.py ${{matrix.options}} -m ${{matrix.momentum}} -p ${{matrix.particle}} # echo "[CI] EVENT DISPLAY SINGLES:" - # event_display d s out/sim.edm4hep.root 0 0 + # event_display d s out/sim.edm4hep.root n 0 0 # mv out/ev out/ev.singles echo "[CI] EVENT DISPLAY INTEGRATED:" - event_display d s out/sim.edm4hep.root + event_display d s out/sim.edm4hep.root n mv out/ev out/ev.all_events echo "[CI] DRAW HITS:" draw_hits d @@ -237,9 +237,70 @@ jobs: retention-days: 1 path: geo*/ + irt_auxfiles: + needs: + - build_3 + runs-on: ubuntu-latest + strategy: + fail-fast: true + matrix: + include: + - { det: d, detector: DRICH } + - { det: p, detector: PFRICH } + steps: + - uses: actions/checkout@v3 + - uses: actions/download-artifact@v3 + with: + name: build_3 + - name: untar_artifacts + run: | + tar xvf build.tar + rm -v build.tar + - uses: cvmfs-contrib/github-action-cvmfs@v3 + - uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "jug_xl:nightly" + setup: environ.sh + run: | + scripts/configure_CI.sh none + echo "[CI] CREATE IRT AUXFILE:" + create_irt_auxfile ${{matrix.det}} geo/${{matrix.detector}}_irt_auxfile.root + - uses: actions/upload-artifact@v3 + with: + name: geometry + retention-days: 1 + path: geo*/ + + sensor_table: + needs: + - build_3 + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: actions/download-artifact@v3 + with: + name: build_3 + - name: untar_artifacts + run: | + tar xvf build.tar + rm -v build.tar + - uses: cvmfs-contrib/github-action-cvmfs@v3 + - uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "jug_xl:nightly" + setup: environ.sh + run: | + scripts/configure_CI.sh none + echo "[CI] MAKE SENSOR TABLE:" + scripts/make_sensor_table.sh + - uses: actions/upload-artifact@v3 + with: + name: geometry + retention-days: 1 + path: geo/sensor_table.* + # EVENT VISUALIZATION ------------------------------------------------------------- -### ### FIXME: not sure how to make these work on CI runners ### - locally: ### run `xvfb-run eic-shell` (to suppress graphics) @@ -316,19 +377,12 @@ jobs: # retention-days: 1 # path: out.*/ -# GENERATE IRT AUXFILES ----------------------------------------------------------- -# - to test it is still possbile +# MISCELLANEOUS TESTS ------------------------------------------------------------- - irt_auxfiles: + test_pixel_gaps: needs: - build_3 runs-on: ubuntu-latest - strategy: - fail-fast: true - matrix: - include: - - { det: d, detector: DRICH } - - { det: p, detector: PFRICH } steps: - uses: actions/checkout@v3 - uses: actions/download-artifact@v3 @@ -345,13 +399,57 @@ jobs: setup: environ.sh run: | scripts/configure_CI.sh none - echo "[CI] CREATE IRT AUXFILE:" - create_irt_auxfile ${{matrix.det}} geo/${{matrix.detector}}_irt_auxfile.root + echo "[CI] RUN SIMULATION:" + simulate.py -s -t4 -n500 + echo "[CI] RUN RECONSTRUCTION:" + recon.rb -c config/recon_digi_only.yaml + echo "[CI] DRAW PIXEL HITS:" + test_pixel_gap_cuts n + - name: tree_artifacts + run: tree out + - name: rename_artifacts + run: mv -v out tests - uses: actions/upload-artifact@v3 with: - name: irt_auxfiles + name: tests retention-days: 1 - path: geo*/ + path: tests*/ + + test_track_propagation: + needs: + - build_3 + runs-on: ubuntu-latest + steps: + - uses: actions/checkout@v3 + - uses: actions/download-artifact@v3 + with: + name: build_3 + - name: untar_artifacts + run: | + tar xvf build.tar + rm -v build.tar + - uses: cvmfs-contrib/github-action-cvmfs@v3 + - uses: eic/run-cvmfs-osg-eic-shell@main + with: + platform-release: "jug_xl:nightly" + setup: environ.sh + run: | + scripts/configure_CI.sh none + echo "[CI] RUN SIMULATION:" + simulate.py -t4 -k10 -n1 + echo "[CI] RUN RECONSTRUCTION:" + recon.rb + echo "[CI] DRAW PIXEL HITS:" + root -b -q scripts/test_tracks.C + - name: tree_artifacts + run: tree out + - name: rename_artifacts + run: mv -v out tests + - uses: actions/upload-artifact@v3 + with: + name: tests + retention-days: 1 + path: tests*/ # ARTIFACT COLLECTION ------------------------------------------------------------- @@ -363,6 +461,9 @@ jobs: # - visual # - optics - irt_auxfiles + - sensor_table + - test_pixel_gaps + - test_track_propagation runs-on: ubuntu-latest strategy: fail-fast: true @@ -382,7 +483,7 @@ jobs: # name: optics - uses: actions/download-artifact@v3 with: - name: irt_auxfiles + name: tests - name: cull run: | find -name .keep | xargs rm -v @@ -393,3 +494,4 @@ jobs: path: | pipeline*/ geo*/ + tests*/ diff --git a/.gitignore b/.gitignore index e99cb061..c04c75ac 100644 --- a/.gitignore +++ b/.gitignore @@ -10,6 +10,7 @@ lib macro !macro/.keep out +geo/*.txt !out/.keep tmp !tmp/.keep diff --git a/Makefile b/Makefile index 0e976fbe..00c1d04e 100644 --- a/Makefile +++ b/Makefile @@ -31,20 +31,25 @@ EICRECON_DIR = ${DRICH_DEV}/EICrecon/src/services/geometry/richgeo DEPS += -I$(EICRECON_DIR) LIB_TARGET = lib -IRTGEO_LIB_NAME = IrtGeo -IRTGEO_LIB = $(LIB_TARGET)/lib$(IRTGEO_LIB_NAME).so -IRTGEO_ROOT = $(EICRECON_DIR) -IRTGEO_SOURCES := $(wildcard $(IRTGEO_ROOT)/IrtGeo*.cc) -IRTGEO_HEADERS := $(wildcard $(IRTGEO_ROOT)/IrtGeo*.h) $(IRTGEO_ROOT)/RichGeo.h +RICHGEO_LIB_NAME = RichGeo +RICHGEO_LIB = $(LIB_TARGET)/lib$(RICHGEO_LIB_NAME).so +RICHGEO_ROOT = $(EICRECON_DIR) +RICHGEO_SOURCES := $(wildcard $(RICHGEO_ROOT)/IrtGeo*.cc) $(RICHGEO_ROOT)/ReadoutGeo.cc +RICHGEO_HEADERS := $(wildcard $(RICHGEO_ROOT)/IrtGeo*.h) $(RICHGEO_ROOT)/ReadoutGeo.h $(RICHGEO_ROOT)/RichGeo.h IRT_AUXFILE_SOURCE = src/create_irt_auxfile.cpp IRT_AUXFILE_EXECUTABLE = $(BIN_TARGET)/create_irt_auxfile SOURCES := $(filter-out $(IRT_AUXFILE_SOURCE), $(SOURCES)) EXECUTABLES := $(filter-out $(IRT_AUXFILE_EXECUTABLE), $(EXECUTABLES)) +PIXEL_GAP_SOURCE = src/test_pixel_gap_cuts.cpp +PIXEL_GAP_EXECUTABLE = $(BIN_TARGET)/test_pixel_gap_cuts +SOURCES := $(filter-out $(PIXEL_GAP_SOURCE), $(SOURCES)) +EXECUTABLES := $(filter-out $(PIXEL_GAP_EXECUTABLE), $(EXECUTABLES)) + #-------------------------------------------- -all: $(EXECUTABLES) $(IRTGEO_LIB) $(IRT_AUXFILE_EXECUTABLE) +all: $(EXECUTABLES) $(RICHGEO_LIB) $(IRT_AUXFILE_EXECUTABLE) $(PIXEL_GAP_EXECUTABLE) $(EXECUTABLES): $(BIN_TARGET)/%: src/%.cpp mkdir -p $(BIN_TARGET) @@ -54,22 +59,34 @@ $(EXECUTABLES): $(BIN_TARGET)/%: src/%.cpp $(CXX) -o $@ $@.o $(LIBS) $(RM) $@.o -$(IRTGEO_LIB): $(IRTGEO_HEADERS) $(IRTGEO_SOURCES) +$(RICHGEO_LIB): $(RICHGEO_HEADERS) $(RICHGEO_SOURCES) ifeq ($(IRT_ROOT_DICT_FOUND),1) mkdir -p $(LIB_TARGET) @echo "----- build $@ -----" - $(CXX) $(IRTGEO_SOURCES) -shared -o $@ $(FLAGS) $(DEPS) $(LIBS) + $(CXX) $(RICHGEO_SOURCES) -shared -o $@ $(FLAGS) $(DEPS) $(LIBS) +else + @echo "WARNING: skip building $@ since IRT ROOT dict not found" +endif + +$(IRT_AUXFILE_EXECUTABLE): $(IRT_AUXFILE_SOURCE) $(RICHGEO_LIB) +ifeq ($(IRT_ROOT_DICT_FOUND),1) + mkdir -p $(BIN_TARGET) + @echo "----- build $@.o -----" + $(CXX) -c $< -o $@.o $(FLAGS) $(DEPS) + @echo "--- make executable $@" + $(CXX) -o $@ $@.o $(LIBS) -L$(LIB_TARGET) -l$(RICHGEO_LIB_NAME) + $(RM) $@.o else @echo "WARNING: skip building $@ since IRT ROOT dict not found" endif -$(IRT_AUXFILE_EXECUTABLE): $(IRT_AUXFILE_SOURCE) $(IRTGEO_LIB) +$(PIXEL_GAP_EXECUTABLE): $(PIXEL_GAP_SOURCE) $(RICHGEO_LIB) ifeq ($(IRT_ROOT_DICT_FOUND),1) mkdir -p $(BIN_TARGET) @echo "----- build $@.o -----" $(CXX) -c $< -o $@.o $(FLAGS) $(DEPS) @echo "--- make executable $@" - $(CXX) -o $@ $@.o $(LIBS) -L$(LIB_TARGET) -l$(IRTGEO_LIB_NAME) + $(CXX) -o $@ $@.o $(LIBS) -L$(LIB_TARGET) -l$(RICHGEO_LIB_NAME) $(RM) $@.o else @echo "WARNING: skip building $@ since IRT ROOT dict not found" @@ -77,5 +94,4 @@ endif clean: @echo "CLEAN ======================================================" - $(RM) $(EXECUTABLES) $(IRT_AUXFILE_EXECUTABLE) $(IRTGEO_LIB) - + $(RM) $(EXECUTABLES) $(IRT_AUXFILE_EXECUTABLE) $(PIXEL_GAP_EXECUTABLE) $(RICHGEO_LIB) diff --git a/README.md b/README.md index 8a13a773..8a8909b5 100644 --- a/README.md +++ b/README.md @@ -392,6 +392,7 @@ flowchart TB - [CERN host](https://root.cern/js/) (recommended) - Local host (advanced, but offers better control) - see [setup guide](doc/jsroot.md) - [ANL hosted](https://eic.phy.anl.gov/geoviewer/) + - alternatively, run `view.py` to use a `ROOT` `TGeoManager` - browse the ROOT file geometry tree in the sidebar on the left: ``` detector_geometry.root diff --git a/config/recon_digi_only.yaml b/config/recon_digi_only.yaml new file mode 100644 index 00000000..aa160c87 --- /dev/null +++ b/config/recon_digi_only.yaml @@ -0,0 +1,33 @@ +################################################################################ +# CONFIGURATION FILE FOR DRICH EICrecon USAGE +# - use with EICrecon `git` branch `irt-algo` or `irt-algo-stable` +################################################################################ + +### PODIO collections to include in the output +podio: + output_include_collections: + ### simulation + - DRICHHits + - MCParticles + ### digitization + - DRICHRawHits + - DRICHRawHitsAssociations + +### DRICH Configuration Overrides +DRICH: + DRICHRawHits: + enablePixelGaps: true + +### EICrecon log levels +log_levels: + eicrecon: info + richgeo: info + DRICH: + DRICHRawHits: info + +### common settings +jana: + nevents: 0 + debug_plugin_loading: 1 +acts: + MaterialMap: calibrations/materials-map.cbor diff --git a/doc/tutorials/1-setup-and-running-simulations.md b/doc/tutorials/1-setup-and-running-simulations.md index 201281a8..df620760 100644 --- a/doc/tutorials/1-setup-and-running-simulations.md +++ b/doc/tutorials/1-setup-and-running-simulations.md @@ -181,11 +181,11 @@ Add the following script, named `jsroot` to your `$PATH`: #!/bin/bash pushd ${HOME}/jsroot # change this to wherever you unpacked the jsroot release echo "FOR FULL DETECTOR:" -printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark\n\n" +printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark&opt=ctrl;\n\n" echo "FOR DRICH ONLY:" -printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark&item=default/world_volume/DRICH_0/DRICH_gas_0&opt=ctrl;\n\n" +printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark&item=default/world_volume/DRICH_0/DRICH_gas_0&opt=ctrl;vislvl8;more\n\n" echo "FOR PFRICH ONLY:" -printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark&item=default/world_volume/PFRICH_0/PFRICH_gas_0&opt=ctrl;\n\n" +printf "\nhttp://localhost:8000/?file=geo/detector_geometry.root&dark&item=default/world_volume/PFRICH_0/PFRICH_gas_0&opt=ctrl;vislvl8;more\n\n" echo "============================================" python -m http.server popd @@ -278,17 +278,13 @@ View the image in `out/ev/`. You may also draw one event at a time (see the usage guide) ```bash -event_display d s out/sim.edm4hep.root 0 0 +event_display d s out/sim.edm4hep.root n 0 0 ``` View the images in `out/ev/`. -Each green box is a single SiPM, and each histogram bin is a single SiPM pixel (8x8 pixels per SiPM). If you want to zoom in, edit `src/event_display.cpp` and uncomment the line -```cpp -#define INTERACTIVE_USE -``` -Then re-compile (`make`) and rerun. This will keep the `TCanvas` open and allow you to zoom in. For example, here is a closeup of the gas ring from a single event, by running +Each green box is a single SiPM, and each histogram bin is a single SiPM pixel (8x8 pixels per SiPM). If you want to zoom in, use interactive mode by changing `n` to `i`, which will keep the `TCanvas` open and allow you to zoom in. For example, here is a closeup of the gas ring from a single event, by running ```bash -event_display d s out/sim.edm4hep.root 0 +event_display d s out/sim.edm4hep.root i 0 ``` ![sim-ev-all-zoom](img/sim-ev-all-zoom.png) diff --git a/scripts/configure_CI.sh b/scripts/configure_CI.sh index 0db37a67..680f866e 100755 --- a/scripts/configure_CI.sh +++ b/scripts/configure_CI.sh @@ -17,7 +17,7 @@ for repo in $repo_list; do git clone https://github.com/eic/irt.git --branch main ;; EICrecon) - git clone https://github.com/eic/EICrecon.git --branch irt-algo + git clone https://github.com/eic/EICrecon.git --branch main ;; reconstruction_benchmarks) git clone https://eicweb.phy.anl.gov/EIC/benchmarks/reconstruction_benchmarks.git --branch irt-algo diff --git a/scripts/draw_sensor_table.py b/scripts/draw_sensor_table.py new file mode 100755 index 00000000..06b2f1ac --- /dev/null +++ b/scripts/draw_sensor_table.py @@ -0,0 +1,34 @@ +#!/usr/bin/env python +# draw sensor positions and orientations from table produced by `make_sensor_table.sh` + +import sys +import ROOT as r + +root_file_name = sys.argv[1] if len(sys.argv)>1 else 'geo/sensor_table.root' +print(f'Opening file {root_file_name}...') +root_file = r.TFile(root_file_name, 'READ') +tr = root_file.Get('table') + +cut = 'sector==0' + +def branch3D(name): + return f'{name}_y:{name}_x:{name}_z' +def cross(a,b): + cx = f'{a}_y*{b}_z-{a}_z*{b}_y' + cy = f'{a}_z*{b}_x-{a}_x*{b}_z' + cz = f'{a}_x*{b}_y-{a}_y*{b}_x' + return f'{cy}:{cx}:{cz}' + +canv = r.TCanvas() +canv.Divide(2,2) +canv.cd(1) +tr.Draw(branch3D('pos'), cut) +canv.cd(2) +tr.Draw(branch3D('normX'), cut) +canv.cd(3) +tr.Draw(branch3D('normY'), cut) +canv.cd(4) +tr.Draw(cross('normX','normY'), cut) # normZ = normX x normY + +root_file.Close() + diff --git a/scripts/make_sensor_table.sh b/scripts/make_sensor_table.sh new file mode 100755 index 00000000..70479f79 --- /dev/null +++ b/scripts/make_sensor_table.sh @@ -0,0 +1,22 @@ +#!/bin/bash +# generate table of sensor positions and orientation vectors + +table_file='geo/sensor_table.txt' +root_file='geo/sensor_table.root' +branches='name/C:ID/C:sector/I:pos_x/D:pos_y/D:pos_z/D:normX_x/D:normX_y/D:normX_z/D:normY_x/D:normY_y/D:normY_z/D' +header=$(echo $branches | sed 's;/.;;g' | sed 's;:; ;g' | sed 's;sector;& ;g' | sed 's;_z;& ;g') +echo "#$header" > $table_file +bin/create_irt_auxfile d | grep sensor_de | sed 's;^.*sensor_de;sensor_de;g' >> $table_file +python -c """ +import ROOT as r +root_file = r.TFile(\"$root_file\", 'RECREATE') +tr = r.TTree('table','table') +tr.ReadFile(\"$table_file\", \"$branches\") +tr.Write() +root_file.Close() +""" +echo """ +CREATED: + $table_file + $root_file +""" diff --git a/scripts/test_pixel_gap_cuts.C b/scripts/test_pixel_gap_cuts.C deleted file mode 100644 index 4770c153..00000000 --- a/scripts/test_pixel_gap_cuts.C +++ /dev/null @@ -1,24 +0,0 @@ -// test if the pixel gap cuts are working -// 1. add member `edm4hep::Vector3d position` to `edm4eic::RawPMTHit` -// 2. fill this member in `PhotoMultiplierHitDigi.cc` -// 3. run a simulation, with varying momentum directions, then run reconstruction -// 4. run this script on the output - -R__LOAD_LIBRARY(fmt) -#include - -void test_pixel_gap_cuts() { - - auto infile = new TFile("out/rec.edm4hep.root"); - auto tr = (TTree*) infile->Get("events"); - - // draw gaps in x,y - if(tr->GetBranch("DRICHRawHits.position.x")!=nullptr) { - auto c = new TCanvas(); - auto h = new TH2D("h","h",10000,-15,15,10000,-15,15); - tr->Project("h","DRICHRawHits.position.y:DRICHRawHits.position.x",""); - fmt::print("ENTRIES = {}\n",h->GetEntries()); - h->Draw(); - c->SaveAs("out/gap.png"); - } -} diff --git a/scripts/test_tracks.C b/scripts/test_tracks.C index de39d821..4f009a8a 100644 --- a/scripts/test_tracks.C +++ b/scripts/test_tracks.C @@ -41,7 +41,7 @@ void test_tracks(TString infileN = "out/rec.edm4hep.root") { double zmin = 1800; double zmax = 3300; - ymax = 200; // zoom to horizontal (y=0) plane + // ymax = 200; // zoom to horizontal (y=0) plane // 2D plots enum views { zx, zy, xy, Nview }; diff --git a/src/event_display.cpp b/src/event_display.cpp index 7b601329..bc1fefa2 100644 --- a/src/event_display.cpp +++ b/src/event_display.cpp @@ -1,12 +1,5 @@ // dRICH event display - -////////////////////////////////// -// if defined, keep TCanvas open for interactive usage (and with extra histograms) -//#define INTERACTIVE_USE -////////////////////////////////// - - // test readout segmentation #include #include @@ -17,17 +10,18 @@ #include // ROOT -#include "TSystem.h" -#include "TStyle.h" -#include "TCanvas.h" -#include "TApplication.h" -#include "TBox.h" -#include "ROOT/RDataFrame.hxx" -#include "ROOT/RDF/HistoModels.hxx" -#include "ROOT/RVec.hxx" +#include +#include +#include +#include +#include +#include +#include +#include // DD4Hep -#include "DD4hep/Detector.h" +#include +#include // local #include "WhichRICH.h" @@ -40,44 +34,42 @@ int main(int argc, char** argv) { // arguments // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ if(argc<=3) { - fmt::print("\nUSAGE: {} [d/p] [s/r] [input_root_file] [event_num_min] [event_num_max]\n\n",argv[0]); + fmt::print("\nUSAGE: {} [d/p] [s/r] [input_root_file] [i/n] [event_num_min] [event_num_max]\n\n",argv[0]); fmt::print(" [d/p]: detector\n"); - fmt::print(" - d for dRICH\n"); - fmt::print(" - p for pfRICH\n"); + fmt::print(" - d: for dRICH\n"); + fmt::print(" - p: for pfRICH\n"); fmt::print("\n"); fmt::print(" [s/r]: file type:\n"); - fmt::print(" - s for simulation file (all photons)\n"); - fmt::print(" - r for reconstructed file (digitized hits)\n"); + fmt::print(" - s: for simulation file (all photons)\n"); + fmt::print(" - r: for reconstructed file (digitized hits)\n"); fmt::print("\n"); fmt::print(" [input_root_file]: output from simulation or reconstruction\n"); fmt::print("\n"); + fmt::print(" [i/n]: interactive mode (optional)\n"); + fmt::print(" - i: keep TCanvases open for interaction\n"); + fmt::print(" - n: do not open TCanvases and just produce images (default)\n"); + fmt::print("\n"); fmt::print(" [event_num_min]: minimum event number (optional)\n"); fmt::print(" - if unspecified, draw sum of all events\n"); fmt::print(" - if specified, but without [event_num_max], draw only\n"); fmt::print(" this event\n"); fmt::print(" - if specified with [event_num_max], draw the range\n"); fmt::print(" of events, ONE at a time\n"); + fmt::print("\n"); fmt::print(" [event_num_max]: maximum event number (optional)\n"); fmt::print(" - set to 0 if you want the maximum possible\n"); fmt::print("\n"); - fmt::print("\n"); - fmt::print("NOTE: INTERACTIVE_USE mode is "); -#ifdef INTERACTIVE_USE - fmt::print("ON: TCanvases (and extra histograms) will remain open\n"); -#else - fmt::print("OFF: TCanvases will saved as PNG files\n"); - fmt::print("- FIXME: there is still a slow memory leak, don't run too many\n"); -#endif - fmt::print("- this setting is hard-coded in {}.cpp\n",argv[0]); return 2; } - std::string zDirectionStr = argv[1]; - std::string fileType = argv[2]; - TString infileN = TString(argv[3]); - int evnumMin = argc>4 ? std::atoi(argv[4]) : -1; - int evnumMax = argc>5 ? std::atoi(argv[5]) : -1; + std::string zDirectionStr = argv[1]; + std::string fileType = argv[2]; + TString infileN = TString(argv[3]); + TString interactiveOpt = argc>4 ? TString(argv[4]) : "n"; + int evnumMin = argc>5 ? std::atoi(argv[5]) : -1; + int evnumMax = argc>6 ? std::atoi(argv[6]) : -1; WhichRICH wr(zDirectionStr); if(!wr.valid) return 1; + bool interactiveOn = interactiveOpt=="i"; // settings // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -92,15 +84,14 @@ int main(int argc, char** argv) { // dilations: for re-scaling module positions and segment positions // for drawing; if you change `numPx`, consider tuning these parameters // as well - const Int_t dilation = 4; + const Double_t dilation = 4.3; // drawing gStyle->SetPalette(55); gStyle->SetOptStat(0); Bool_t singleCanvas = true; // if true, draw all hitmaps on one canvas -#ifndef INTERACTIVE_USE - singleCanvas = true; -#endif + if(!interactiveOn) + singleCanvas = true; // data collections std::string inputCollection; @@ -117,10 +108,10 @@ int main(int argc, char** argv) { // setup // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ -#ifdef INTERACTIVE_USE - // define application environment, to keep canvases open - TApplication mainApp("mainApp",&argc,argv); -#endif + // define application environment, to keep canvases open, if interactive mode is on + TApplication *mainApp; + if(interactiveOn) + mainApp = new TApplication("mainApp",&argc,argv); // main dataframe RDataFrame dfIn("events",infileN.Data()); @@ -146,10 +137,10 @@ int main(int argc, char** argv) { fmt::print("Reading event numbers {} to {}, one at a time\n",evnumMin,evnumMax); for(int e=evnumMin; e<=evnumMax; e++) evnumRanges.push_back({ e, e+1 }); -#ifdef INTERACTIVE_USE - fmt::print(stderr,"ERROR: cannot yet run with INTERACTIVE_USE on an event number range; change it in {}.cpp and rebuild\n",argv[0]); - return 1; -#endif + if(interactiveOn) { + fmt::print(stderr,"ERROR: cannot yet run with interactive mode enabled with an event number range (FIXME)\n"); + return 1; + } } @@ -167,8 +158,6 @@ int main(int argc, char** argv) { const auto det = &(Detector::getInstance()); det->fromXML(compactFile); const auto detRich = det->detector(richName); - const auto posRich = detRich.placement().position(); - const auto cellMask = ULong_t(std::stoull(det->constant(wr.XRICH+"_cell_mask"))); const auto nSectors = wr.zDirection>0 ? det->constant(wr.XRICH+"_num_sectors") : 1; // cellID decoder @@ -196,29 +185,30 @@ int main(int argc, char** argv) { }; // build sensor position LUT `imod2hitmapXY` - /* - find the sector 0 sensors, and build a map of their module number `imod` to + /* - find the sector 0 sensors, and build a map of their module number `pdu`,`sipm` to * X and Y coordinates to use in the hitmap * - these hitmap coordinates are from the sensor position X and Y, rescaled * by the factor `dilation` and rounded to the nearest integer * - also builds a list of `TBox`es, for drawing the sensors on the hitmap - * - the unique ID of the sensor `Detector`, called `imodsec`, includes `imod` - * - `cellID & cellMask` should be equivalent to `imodsec`; therefore, - * `imodsec` can be converted to `imod` by decoding `imodsec` the same way - * we would decode `cellID` */ - std::map> imod2hitmapXY; + std::map,std::pair> imod2hitmapXY; std::vector boxList; for(auto const& [de_name, detSensor] : detRich.children()) { if(de_name.find(wr.sensorNamePattern)!=std::string::npos) { // convert global position to hitmapX and Y - auto posSensor = posRich + detSensor.placement().position(); - auto hitmapX = Long64_t(dilation*posSensor.x() + 0.5); - auto hitmapY = Long64_t(dilation*posSensor.y() + 0.5); + const auto detSensorPars = detSensor.extension(true); + if(detSensorPars==nullptr) + throw std::runtime_error(fmt::format("sensor '{}' does not have VariantParameters", de_name)); + auto posSensorX = detSensorPars->get("pos_x"); + auto posSensorY = detSensorPars->get("pos_y"); + auto hitmapX = Long64_t(dilation*posSensorX + 0.5); + auto hitmapY = Long64_t(dilation*posSensorY + 0.5); // convert unique cellID to module number, using the cellID decoder - auto imodsec = ULong_t(detSensor.id()); - auto imod = decodeCellID("module")(RVecUL({imodsec})).front(); + auto sensorID = ULong_t(detSensor.id()); + auto pdu = decodeCellID("pdu")(RVecUL({sensorID})).front(); + auto sipm = decodeCellID("sipm")(RVecUL({sensorID})).front(); // add to `imod2hitmapXY` and create sensor `TBox` - imod2hitmapXY.insert({imod,{hitmapX,hitmapY}}); + imod2hitmapXY.insert({{pdu,sipm}, {hitmapX,hitmapY}}); boxList.push_back(new TBox( hitmapX, hitmapY, @@ -231,26 +221,30 @@ int main(int argc, char** argv) { } } - // convert vector of `imod`s to vector of hitmap X or Y - auto imod2hitmapXY_get = [&imod2hitmapXY] (RVecL imodVec, int c) { + // convert vectors of `pdu`, `sipm` to vector of hitmap X or Y + auto imod2hitmapXY_get = [&imod2hitmapXY] (RVecL pduVec, RVecL sipmVec, int c) { RVecL result; Long64_t pos; - for(auto imod : imodVec) { + if(pduVec.size() != sipmVec.size()) + throw std::runtime_error("imod2hitmapXY_get input vectors differ in size"); + for(unsigned long i=0; iDivide(3,2); - int pad=1; - for(auto hist : fieldHists) { - c->GetPad(pad)->SetLogy(); + if(interactiveOn) { + // draw cellID field histograms + c = new TCanvas(); + c->Divide(4,2); + int pad=1; + for(auto hist : fieldHists) { + c->GetPad(pad)->SetLogy(); + c->cd(pad); + // if(TString(hist->GetName())!="pdu") hist->SetBarWidth(4); + hist->SetLineColor(kBlack); + hist->SetFillColor(kBlack); + hist->Draw("bar"); + pad++; + } + + // draw segmentation XY plot, along with expected box c->cd(pad); - if(TString(hist->GetName())!="module") hist->SetBarWidth(4); - hist->SetLineColor(kBlack); - hist->SetFillColor(kBlack); - hist->Draw("bar"); - pad++; + c->GetPad(pad)->SetGrid(1,1); + segXY->Draw("colz"); + auto expectedBox = new TBox(segXmin, segXmin, segXmax+1, segXmax+1); + expectedBox->SetFillStyle(0); + expectedBox->SetLineColor(kBlack); + expectedBox->SetLineWidth(4); + expectedBox->Draw("same"); } - // draw segmentation XY plot, along with expected box - c->cd(pad); - c->GetPad(pad)->SetGrid(1,1); - segXY->Draw("colz"); - auto expectedBox = new TBox(segXmin, segXmin, segXmax+1, segXmax+1); - expectedBox->SetFillStyle(0); - expectedBox->SetLineColor(kBlack); - expectedBox->SetLineWidth(4); - expectedBox->Draw("same"); -#endif - // draw pixel hitmap if(singleCanvas) { c = new TCanvas("canv","canv",3*700,2*600); c->Divide(3,2); }; int secBin; @@ -394,24 +390,23 @@ int main(int argc, char** argv) { }; // either hold the TCanvases open, or save them as PNG files -#ifdef INTERACTIVE_USE - fmt::print("\n\npress ^C to exit.\n\n"); - mainApp.Run(); -#else - gROOT->ProcessLine(".! mkdir -p out/ev"); - TString pngName = Form("out/ev/%s.png", fmt::format("{:08}",evnum).c_str()); - c->SaveAs(pngName); - fmt::print("Saved image: {}", pngName); - // cleanup and avoid memory leaks # FIXME: refactor this... and there is still a slow leak... - delete c; - for(int sec=0; secRun(); + } else { + gROOT->ProcessLine(".! mkdir -p out/ev"); + TString pngName = Form("out/ev/%s.png", fmt::format("{:08}",evnum).c_str()); + c->SaveAs(pngName); + fmt::print("Saved image: {}", pngName); + // cleanup and avoid memory leaks # FIXME: refactor this... and there is still a slow leak... + delete c; + for(int sec=0; sec + +#include +#include +#include +#include + +#include +#include + +#include + +#include + +int main(int argc, char** argv) { + + // args + if(argc<=1) { + fmt::print("USAGE: {} [n/i] [file(default=out/rec.edm4hep.root)]\n", argv[0]); + fmt::print(" [n/i]: n=non-interactive, i=interactive\n"); + return 2; + } + std::string interactiveOpt = std::string(argv[1]); + std::string root_file_name = argc>2 ? std::string(argv[2]) : "out/rec.edm4hep.root"; + + // set interactive mode + bool interactiveOn = interactiveOpt=="i"; + TApplication *app; + if(interactiveOn) + app = new TApplication("app", &argc, argv); + + // logger + auto logger = spdlog::default_logger()->clone("richgeo"); + logger->set_level(spdlog::level::trace); + + // geometry from compact file + std::string DETECTOR_PATH(getenv("DETECTOR_PATH")); + std::string DETECTOR_CONFIG(getenv("DETECTOR_CONFIG")); + if(DETECTOR_PATH.empty() || DETECTOR_CONFIG.empty()) { + logger->error("cannot find default compact file, since env vars DETECTOR_PATH and DETECTOR_CONFIG are not set"); + return 1; + } + auto compactFile = DETECTOR_PATH + "/" + DETECTOR_CONFIG + ".xml"; + auto det = &dd4hep::Detector::getInstance(); + det->fromXML(compactFile); + + // ReadoutGeo + richgeo::ReadoutGeo geo("DRICH", det, logger); + + // open input file + auto reader = podio::ROOTFrameReader(); + reader.openFile(root_file_name); + const std::string tree_name = "events"; + + // local hits histogram + auto h = new TH2D("h","local MC SiPM hits",10000,-15,15,10000,-15,15); + + // event loop + for(unsigned e=0; etrace("EVENT {}", e); + auto frame = podio::Frame(reader.readNextEntry(tree_name)); + + const auto& hit_assocs = frame.get("DRICHRawHitsAssociations"); + if(!hit_assocs.isValid()) + throw std::runtime_error("cannot find hit associations"); + + for(const auto& hit_assoc : hit_assocs) { + for(const auto& sim_hit : hit_assoc.getSimHits()) { + + auto cellID = sim_hit.getCellID(); + auto pos = sim_hit.getPosition(); + dd4hep::Position pos_global(pos.x*dd4hep::mm, pos.y*dd4hep::mm, pos.z*dd4hep::mm); + auto pos_local = geo.GetSensorLocalPosition(cellID, pos_global); + h->Fill(pos_local.y()/dd4hep::mm, pos_local.x()/dd4hep::mm); + + if(logger->level() <= spdlog::level::trace) { + logger->trace("pixel hit on cellID={:#018x}",cellID); + auto print_pos = [&] (std::string name, dd4hep::Position p) { + logger->trace(" {:>30} x={:.2f} y={:.2f} z={:.2f} [mm]: ", name, p.x()/dd4hep::mm, p.y()/dd4hep::mm, p.z()/dd4hep::mm); + }; + print_pos("pos_local", pos_local); + print_pos("pos_global", pos_global); + // auto dim = m_cellid_converter->cellDimensions(cellID); + // for (std::size_t j = 0; j < std::size(dim); ++j) + // logger->trace(" - dimension {:<5} size: {:.2}", j, dim[j]); + } + } + } + } + + gStyle->SetOptStat(0); + auto c = new TCanvas(); + h->Draw(); + fmt::print("NUMBER OF DIGITIZED PHOTONS: {}\n", h->GetEntries()); + if(interactiveOn) { + fmt::print("\n\npress ^C to exit.\n\n"); + app->Run(); + } else { + c->SaveAs("out/pixel_gaps.png"); + } +} diff --git a/view.py b/view.py new file mode 100755 index 00000000..e7ccc852 --- /dev/null +++ b/view.py @@ -0,0 +1,8 @@ +#!/usr/bin/env python +import ROOT as r +r.TGeoManager.Import("geo/detector_geometry.root") +r.gGeoManager.SetVisLevel(8) +drich = r.gGeoManager.GetVolume("DRICH") +drich.Draw("ogl") +# sensor = r.gGeoManager.GetVolume("DRICH_sensor_sec0") +# sensor.Draw("ogl")