diff --git a/.github/workflows/validation.yml b/.github/workflows/validation.yml index e1f143143..40bd20705 100644 --- a/.github/workflows/validation.yml +++ b/.github/workflows/validation.yml @@ -213,72 +213,15 @@ jobs: with: key: ${{ env.BRANCH_NAME }}-${{ github.sha }} path: ${{ env.REST_PATH }} - - name: 02_PandaXiiiMC + - name: PandaXIII Geant4 run: | source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC restG4 ${{ env.REST_PATH }}/examples/restG4/05.PandaXIII/Xe136bb0n.rml -o Xe136bb0n.root - restRoot -b -q ${{ env.REST_PATH }}/examples/restG4/05.PandaXIII/Validate.C'("Xe136bb0n.root")' - # - name: 02_PandaXiiiMC small for reference # We can use these artifacts to update the validation reference file - # run: | - # source ${{ env.REST_PATH }}/thisREST.sh - # cd framework/pipeline/pandaxiii_MC - # restG4 ${{ env.REST_PATH }}/examples/restG4/05.PandaXIII/Xe136bb0n.rml -o Xe136bb0n_small_reference_pipeline.root -n 5 - # restManager --c AllProcesses.rml --i Xe136bb0n_small_reference_pipeline.root --o Xe136bb0n_small_reference_processed_pipeline.root - # restRoot -b -q ../MakeBasicTree.C'("Xe136bb0n_small_reference_processed_pipeline.root", "Xe136bb0n_validation_pipeline.root")' - # - name: Upload Artifacts - # uses: actions/upload-artifact@v3 - # with: - # name: Xe136bb - # path: | - # framework/pipeline/pandaxiii_MC/Xe136bb0n.root - # framework/pipeline/pandaxiii_MC/Xe136bb0n_small_reference_pipeline.root - # framework/pipeline/pandaxiii_MC/Xe136bb0n_small_reference_processed_pipeline.root - # framework/pipeline/pandaxiii_MC/Xe136bb0n_validation_pipeline.root - # retention-days: 1 - - name: PandaXIII Topological - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC - restManager --c AllProcesses.rml --i Xe136bb0n_small_reference.root --o Xe136bb0n_small_reference_processed.root - restRoot -b -q ../MakeBasicTree.C'("Xe136bb0n_small_reference_processed.root")' - restRoot -b -q ../ValidateTrees.C'("Xe136bb0n_validation.root")' - # restManager --c plots.rml --i Xe136bb0n_small_reference_processed.root - # This command is failing in the docker, reproducible locally, not clear why... - # echo | sleep 5 | restManager --c plots.rml --i testOutput.root - # - name: Upload Artifacts - # uses: actions/upload-artifact@v3 - # with: - # name: PandaTrackParam - # path: framework/pipeline/pandaxiii_MC/trackParameter.png - # retention-days: 1 - - name: PandaXIII Topological from Geant4 - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC - echo "using just-generated g4 file" - restManager --c AllProcesses.rml --i Xe136bb0n.root --o Xe136bb0n_processed.root --j 1 --e 5 - restRoot -b -q ../MakeBasicTree.C'("Xe136bb0n_processed.root")' - restRoot -b -q ../ValidateTrees.C'("Xe136bb0n_validation.root")' - - name: Upload Artifacts - uses: actions/upload-artifact@v3 - with: - name: PandaTestOutput - path: framework/pipeline/pandaxiii_MC/Xe136bb0n_processed.root - retention-days: 1 - name: PandaXIII Data run: | source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_data + cd framework/pipeline/panda-x restManager --c P3AutoChain.rml --i CoBo_AsAd0_2019-03-15.graw --o testOutput.root --j 1 - restRoot -b -q ../MakeBasicTree.C'("testOutput.root")' - restRoot -b -q ../ValidateTrees.C'("validation.root")' - - name: Upload Artifacts - uses: actions/upload-artifact@v3 - with: - name: PandaTriggerRate - path: framework/pipeline/pandaxiii_data/TriggerRate.png - retention-days: 1 trex-dm: name: "TREX-DM" @@ -455,77 +398,6 @@ jobs: key: ${{ env.BRANCH_NAME }}-${{ github.sha }} path: ${{ env.REST_PATH }} - - pandax-iii-reference: - name: "PandaX-III on reference version" - runs-on: ubuntu-latest - container: - image: ghcr.io/lobis/root-geant4-garfield:rest-for-physics-reference-jun2022 - needs: [ framework-install-reference ] - steps: - - uses: actions/checkout@v3 - with: - repository: rest-for-physics/framework - path: framework - - name: Checkout framework branch - run: | - cd framework - ./scripts/checkoutRemoteBranch.sh ${{ env.BRANCH_NAME }} - - name: Restore cache - uses: actions/cache@v3 - id: framework-install-cache-reference - with: - key: ${{ env.BRANCH_NAME }}-${{ github.sha }} - path: ${{ env.REST_PATH }} - - name: 02_PandaXiiiMC - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC - restG4 ${{ env.REST_PATH }}/examples/restG4/05.PandaXIII/Xe136bb0n.rml -o Xe136bb0n.root - restRoot -b -q ${{ env.REST_PATH }}/examples/restG4/05.PandaXIII/Validate.C'("Xe136bb0n.root")' - - name: Upload Artifacts - uses: actions/upload-artifact@v3 - with: - name: Xe136bbRef - path: framework/pipeline/pandaxiii_MC/Xe136bb0n.root - retention-days: 1 - - name: PandaXIII Topological - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC - restManager --c AllProcesses.rml --i Xe136bb0n_small_reference.root --o Xe136bb0n_small_reference_processed.root --j 1 --e 5 - restRoot -b -q ../MakeBasicTree.C'("Xe136bb0n_small_reference_processed.root")' - restRoot -b -q ../ValidateTrees.C'("Xe136bb0n_validation.root")' - # restManager --c plots.rml --i Xe136bb0n_small_reference_processed.root - - name: PandaXIII Topological from Geant4 - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_MC - echo "using just-generated g4 file" - restManager --c AllProcesses.rml --i Xe136bb0n.root --o Xe136bb0n_processed.root --j 1 --e 5 - restRoot -b -q ../MakeBasicTree.C'("Xe136bb0n_processed.root")' - # TODO: regenerate reference file - # restRoot -b -q ../ValidateTrees.C'("Xe136bb0n_validation_G4_v10.4.3.root")' - - name: Upload Artifacts - uses: actions/upload-artifact@v3 - with: - name: PandaTestOutputRef - path: framework/pipeline/pandaxiii_MC/Xe136bb0n_processed.root - retention-days: 1 - - name: PandaXIII Data - run: | - source ${{ env.REST_PATH }}/thisREST.sh - cd framework/pipeline/pandaxiii_data - restManager --c P3AutoChain.rml --i CoBo_AsAd0_2019-03-15.graw --o testOutput.root --j 1 - restRoot -b -q ../MakeBasicTree.C'("testOutput.root")' - restRoot -b -q ../ValidateTrees.C'("validation.root")' - - name: Upload Artifacts - uses: actions/upload-artifact@v3 - with: - name: PandaTriggerRateRef - path: framework/pipeline/pandaxiii_data/TriggerRate.png - retention-days: 1 - examples-reference: name: Run examples on reference version runs-on: ubuntu-latest diff --git a/cmake/thisREST.cmake b/cmake/thisREST.cmake index 54566035d..d0cd41b62 100644 --- a/cmake/thisREST.cmake +++ b/cmake/thisREST.cmake @@ -16,6 +16,8 @@ execute_process( string(REGEX REPLACE "\n$" "" GEANT4_PATH "${GEANT4_PATH}") get_filename_component(GEANT4_BIN_DIR "${GEANT4_PATH}/bin/" REALPATH) +set(g4LibPath "") +set(loadG4 "") if (${REST_G4} MATCHES "ON") # https://github.com/rest-for-physics/framework/issues/331 set(g4LibPath ":${GEANT4_PATH}/lib/") @@ -23,15 +25,16 @@ if (${REST_G4} MATCHES "ON") "\# if geant4.sh script is found we load the same Geant4 version as used in compilation\nif [[ -f \\\"${GEANT4_BIN_DIR}/geant4.sh\\\" ]]; then [[ -n \\\"\\\${ZSH_VERSION}\\\" ]] && pushd ${GEANT4_BIN_DIR} > /dev/null\n source ${GEANT4_BIN_DIR}/geant4.sh\n [[ -n \\\"\\\${ZSH_VERSION}\\\" ]] && popd > /dev/null\nfi\n" ) -else () - set(g4LibPath "") - set(loadG4 "") -endif (${REST_G4} MATCHES "ON") +endif () +set(loadMPFR "") if (DEFINED MPFR_PATH) set(loadMPFR "export LD_LIBRARY_PATH=${MPFR_PATH}/lib:\$LD_LIBRARY_PATH") -else () - set(loadMPFR "") +endif () + +set(loadCRY "") +if (DEFINED REST_CRY_PATH) + set(loadCRY "export LD_LIBRARY_PATH=${REST_CRY_PATH}/lib:\$LD_LIBRARY_PATH") endif () set(loadGarfield "") @@ -111,6 +114,7 @@ fi ${loadG4} ${loadMPFR} +${loadCRY} ${loadGarfield} if [ \\\$REST_PATH ] ; then @@ -124,6 +128,8 @@ fi export REST_SOURCE=${CMAKE_CURRENT_SOURCE_DIR} export REST_PATH=\\\${thisdir} +# REST_HOME is where temporary files are stored +export REST_HOME=\\\${HOME} export ROOT_INCLUDE_PATH=\\\$REST_PATH/include${Garfield_INCLUDE_ENV}:\\\$ROOT_INCLUDE_PATH export REST_INPUTDATA=\\\$REST_PATH/data export REST_GARFIELD_INCLUDE=${Garfield_INCLUDE_DIRS} @@ -137,6 +143,7 @@ export PYTHONPATH=${PYTHON_BINDINGS_INSTALL_DIR}:\\\$PYTHONPATH alias restRoot=\\\"restRoot -l\\\" alias restRootMacros=\\\"restRoot -l --m\\\" +alias restListMacros=\\\"restManager ListMacros\\\" if [ \\\$(rest-config --flags | grep \\\"REST_WELCOME=ON\\\") ]; then rest-config --welcome diff --git a/doc/doxygen/images/drawGainMap.png b/doc/doxygen/images/drawGainMap.png new file mode 100644 index 000000000..f6d3df52a Binary files /dev/null and b/doc/doxygen/images/drawGainMap.png differ diff --git a/doc/doxygen/images/drawSpectrum.png b/doc/doxygen/images/drawSpectrum.png new file mode 100644 index 000000000..3d5e148e7 Binary files /dev/null and b/doc/doxygen/images/drawSpectrum.png differ diff --git a/doc/doxygen/images/gainCorrectionComparison.png b/doc/doxygen/images/gainCorrectionComparison.png new file mode 100644 index 000000000..d1dea811c Binary files /dev/null and b/doc/doxygen/images/gainCorrectionComparison.png differ diff --git a/examples/01.alphaTrack/README.md b/examples/01.alphaTrack/README.md index 823f6345b..814af17bf 100644 --- a/examples/01.alphaTrack/README.md +++ b/examples/01.alphaTrack/README.md @@ -22,7 +22,7 @@ This will produce an initial ROOT file that contains the simulated Geant4 event In order to further process the data we need to execute: ``` -restManager --c processing --f data/Run_5MeV_1um.root +restManager --c processing.rml --f data/Run_5MeV_1um.root ``` The processing includes the electron diffusion, the readout segmentation, and shaping and noise signal effects. diff --git a/examples/calibrationCorrection.rml b/examples/calibrationCorrection.rml new file mode 100644 index 000000000..c233b5bd9 --- /dev/null +++ b/examples/calibrationCorrection.rml @@ -0,0 +1,24 @@ + + + + + + + + + + + + + + + + diff --git a/macros/REST_CheckRunFileList.C b/macros/REST_CheckValidRuns.C similarity index 55% rename from macros/REST_CheckRunFileList.C rename to macros/REST_CheckValidRuns.C index eba5324da..2031a7eb9 100644 --- a/macros/REST_CheckRunFileList.C +++ b/macros/REST_CheckValidRuns.C @@ -1,3 +1,4 @@ +#include #include #include @@ -5,20 +6,41 @@ #include "TRestTask.h" #include "TSystem.h" using namespace std; +namespace fs = std::filesystem; -#ifndef RESTTask_CheckRunFileList -#define RESTTask_CheckRunFileList +#ifndef RESTTask_CheckValidRuns +#define RESTTask_CheckValidRuns //******************************************************************************************************* //*** -//*** Your HELP is needed to verify, validate and document this macro -//*** This macro might need update/revision. +//*** Description: This macro will identify run files that were not properly closed. +//*** +//*** -------------- +//*** The first and mandatory argument must provide a full data and pattern to filter the files that +//*** will be checked. E.g. to add all the Hits files from run number 123 the first argument would be +//*** pattern = "/path/to/data/Run00123*Hits*root". +//*** +//*** IMPORTANT: The pattern must be given using double quotes "" +//*** +//*** -------------- +//*** Usage: restManager CheckValidRuns "/full/path/file_*pattern*.root" [purge] +//*** -------------- +//*** +//*** An optional parameter `purge` can be made true. In that case, this macro will not just provide +//*** a list of the files not properly closed but it will also remove them! +//*** +//*** The following command will remove the non-valid runs +//*** -------------- +//*** Usage: restManager CheckValidRuns "/full/path/file_*pattern*.root" 1 +//*** -------------- +//*** +//*** CAUTION: Be aware that any non-REST file in the list will be removed if you use purge=1 //*** //******************************************************************************************************* -Int_t REST_CheckRunFileList(TString namePattern, Int_t N = 100000) { +Int_t REST_CheckValidRuns(TString namePattern, Bool_t purge = false) { TGeoManager::SetVerboseLevel(0); - vector filesNotWellClosed; + vector filesNotWellClosed; TRestStringOutput RESTLog; @@ -70,7 +92,16 @@ Int_t REST_CheckRunFileList(TString namePattern, Int_t N = 100000) { RESTLog << "---------------------" << RESTendl; RESTLog << "Files not well closed" << RESTendl; RESTLog << "---------------------" << RESTendl; - for (int i = 0; i < filesNotWellClosed.size(); i++) RESTLog << filesNotWellClosed[i] << RESTendl; + for (int i = 0; i < filesNotWellClosed.size(); i++) { + RESTLog << filesNotWellClosed[i] << RESTendl; + if (purge) fs::remove(filesNotWellClosed[i]); + } + + if (purge) { + RESTLog << "---------------------------" << RESTendl; + RESTLog << "The files have been removed" << RESTendl; + RESTLog << "---------------------------" << RESTendl; + } } RESTLog << "------------------------------" << RESTendl; diff --git a/macros/REST_GenerateDataSets.C b/macros/REST_GenerateDataSets.C index d6c3f02bc..3c468ca6d 100644 --- a/macros/REST_GenerateDataSets.C +++ b/macros/REST_GenerateDataSets.C @@ -15,14 +15,16 @@ //*** //******************************************************************************************************* -Int_t REST_GenerateDataSets(const std::string& inputRML, const std::string& datasets) { +Int_t REST_GenerateDataSets(const std::string& inputRML, const std::string& datasets, + std::string outPath = "") { std::vector sets = REST_StringHelper::Split(datasets, ","); for (const auto& set : sets) { std::cout << "Set : " << set << std::endl; TRestDataSet d(inputRML.c_str(), set.c_str()); d.GenerateDataSet(); - d.Export("Dataset_" + set + ".root"); + if (!outPath.empty() && outPath.back() != '/') outPath += '/'; + d.Export(outPath + "Dataset_" + set + ".root"); } return 0; } diff --git a/macros/REST_OpenInputFile.C b/macros/REST_OpenInputFile.C index 1bf55e0a5..3087791c2 100644 --- a/macros/REST_OpenInputFile.C +++ b/macros/REST_OpenInputFile.C @@ -46,8 +46,10 @@ void REST_OpenInputFile(const std::string& fileName) { printf("\n%s\n", " - dSet->PrintMetadata()"); printf("%s\n", " - dSet->GetDataFrame().GetColumnNames()"); printf("%s\n\n", " - dSet->GetTree()->GetEntries()"); + printf("%s\n\n", " - dSet->GetDataFrame().Display({\"colName1,colName2\"})->Print()"); if (dSet) delete dSet; dSet = new TRestDataSet(); + dSet->EnableMultiThreading(false); dSet->Import(fileName); } else { printf("\n%s is not a valid TRestRun or TRestDataSet\n", fileName.c_str()); diff --git a/pipeline/pandaxiii_data/CoBo_AsAd0_2019-03-15.graw b/pipeline/panda-x/CoBo_AsAd0_2019-03-15.graw similarity index 100% rename from pipeline/pandaxiii_data/CoBo_AsAd0_2019-03-15.graw rename to pipeline/panda-x/CoBo_AsAd0_2019-03-15.graw diff --git a/pipeline/panda-x/P3AutoChain.rml b/pipeline/panda-x/P3AutoChain.rml new file mode 100644 index 000000000..75327d268 --- /dev/null +++ b/pipeline/panda-x/P3AutoChain.rml @@ -0,0 +1,91 @@ + + + + + + + + //change this value to "auto" to enable database + + + + + + + + + + + + + + + /// We define all observables except MinValue because it is not yet in validation chain + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/pipeline/pandaxiii_MC/AllProcesses.rml b/pipeline/pandaxiii_MC/AllProcesses.rml deleted file mode 100644 index 04a7329bc..000000000 --- a/pipeline/pandaxiii_MC/AllProcesses.rml +++ /dev/null @@ -1,68 +0,0 @@ - - - - - - - - - %options are : silent, warning, info, debug - - // set to 0 as random seed - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/pandaxiii_MC/README.md b/pipeline/pandaxiii_MC/README.md deleted file mode 100644 index e33796ef4..000000000 --- a/pipeline/pandaxiii_MC/README.md +++ /dev/null @@ -1,13 +0,0 @@ -The following scripts are used to validate the PandaX-III simulation & topological analysis chain. - -- **step 1**: Use `restG4 Xe136bb0n.rml` to run test Geant4 simulation. Use `ValidateG4.C` to check if the simulation - produced the expected output. i.e. all 500 events are simulated and metadata is correctly recorded. The output file - will be used later on. - -- **step 2**: Use `AllProcesses.rml` to run restManager for REST simulation and topological analysis. `processes_2D.rml` - and `readout-140kg.rml` is included by it. - -- **step 3**: Use `plots.rml` to check the plots can be successfully drawn by TRestAnalysisPlot - -File `Xe136bb0n_small_reference.root` corresponds to the file generated by step 1 using Geant4 11.0.2 but with a smaller -number of simulated events to reduce size. events. diff --git a/pipeline/pandaxiii_MC/Xe136bb0n_small_reference.root b/pipeline/pandaxiii_MC/Xe136bb0n_small_reference.root deleted file mode 100644 index ea364d450..000000000 Binary files a/pipeline/pandaxiii_MC/Xe136bb0n_small_reference.root and /dev/null differ diff --git a/pipeline/pandaxiii_MC/Xe136bb0n_validation.root b/pipeline/pandaxiii_MC/Xe136bb0n_validation.root deleted file mode 100644 index 300359553..000000000 Binary files a/pipeline/pandaxiii_MC/Xe136bb0n_validation.root and /dev/null differ diff --git a/pipeline/pandaxiii_MC/Xe136bb0n_validation_G4_v10.4.3.root b/pipeline/pandaxiii_MC/Xe136bb0n_validation_G4_v10.4.3.root deleted file mode 100644 index 4deb656af..000000000 Binary files a/pipeline/pandaxiii_MC/Xe136bb0n_validation_G4_v10.4.3.root and /dev/null differ diff --git a/pipeline/pandaxiii_MC/Xe136bb0n_validation_root6.20.root b/pipeline/pandaxiii_MC/Xe136bb0n_validation_root6.20.root deleted file mode 100644 index e63912f3e..000000000 Binary files a/pipeline/pandaxiii_MC/Xe136bb0n_validation_root6.20.root and /dev/null differ diff --git a/pipeline/pandaxiii_MC/plots.rml b/pipeline/pandaxiii_MC/plots.rml deleted file mode 100644 index efc5ec537..000000000 --- a/pipeline/pandaxiii_MC/plots.rml +++ /dev/null @@ -1,74 +0,0 @@ - - - - - - //change this value to "auto" to enable database - //change this value to "auto" to enable database - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/pandaxiii_MC/processes_2D.rml b/pipeline/pandaxiii_MC/processes_2D.rml deleted file mode 100644 index aa6aa3cb0..000000000 --- a/pipeline/pandaxiii_MC/processes_2D.rml +++ /dev/null @@ -1,587 +0,0 @@ - - - - - - - - - - - - - - - - // fwhm - // absolute gain // parameter name="electronicsGain" - value="671744" // electrons for 4096 ADC units - - - - - - - - - // fwhm - // absolute gain // parameter name="electronicsGain" - value="671744" // electrons for 4096 ADC units - - - - - - - - - - - - - - - - - - // Number of sigmas to perform the calculation - - - - // // electrons in each time bin - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - //response file to be used - to shape the signal - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/pandaxiii_MC/readout-140kg.rml b/pipeline/pandaxiii_MC/readout-140kg.rml deleted file mode 100644 index e4ff8fc72..000000000 --- a/pipeline/pandaxiii_MC/readout-140kg.rml +++ /dev/null @@ -1,51 +0,0 @@ - - - - //The Microbulk MicroMegas - - - // Y-strips (pixel has fixed x coordinate) - - - // Last strip is special - - - // X-strips (pixel has fixed y coordinate) - - // First strip is special - - - - - - - - - - // 4 modules - - - // 6-modules - - - // 8-modules - - - // 8-modules - - - // 8-modules - - - // 8-modules - - - // 6-modules - - - // 4-modules - - - - - diff --git a/pipeline/pandaxiii_MC/readouts.rml b/pipeline/pandaxiii_MC/readouts.rml deleted file mode 100644 index 09838438b..000000000 --- a/pipeline/pandaxiii_MC/readouts.rml +++ /dev/null @@ -1,342 +0,0 @@ - - - - - - - - - - - - - - - //The Microbulk MicroMegas - - - // Y-strips (pixel has fixed x coordinate) - - - - - - - - - - // Last strip is special - - - - - - - - - // X-strips (pixel has fixed y coordinate) - - // First strip is special - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - // 4 modules - - - - - // 6-modules - - - - - // 7-modules - - - - - // 7-modules - - - - - // 7-modules - - - - - // 6-modules - - - - - // 4-modules - - - - - - - - - - - - - - - - - - - - - - - // 2 modules - - - - - // 5-modules - - - - - // 6-modules - - - - - // 7-modules - - - - - // 6-modules - - - - - // 5-modules - - - - - // 2-modules - - - - - - - - diff --git a/pipeline/pandaxiii_data/P3AutoChain.rml b/pipeline/pandaxiii_data/P3AutoChain.rml deleted file mode 100644 index 637f2d5b7..000000000 --- a/pipeline/pandaxiii_data/P3AutoChain.rml +++ /dev/null @@ -1,90 +0,0 @@ - - - - - - - - //change this value to "auto" to enable database - - - - - - - - - - - - - - - /// We define all observables except MinValue because it is not yet in validation chain - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - diff --git a/pipeline/pandaxiii_data/validation.root b/pipeline/pandaxiii_data/validation.root deleted file mode 100644 index 5d74aa450..000000000 Binary files a/pipeline/pandaxiii_data/validation.root and /dev/null differ diff --git a/pipeline/validateMacros.py b/pipeline/validateMacros.py index 0b3d249ae..7a7520320 100755 --- a/pipeline/validateMacros.py +++ b/pipeline/validateMacros.py @@ -1,6 +1,6 @@ import os, sys -os.system("restRoot -b -q -m 1 > error.log 2>&1") +os.system("restRoot -b -q --m 1 > error.log 2>&1") with open("error.log") as f: lines = f.readlines() diff --git a/source/bin/restRoot.cxx b/source/bin/restRoot.cxx index e59bfb97d..17e2a1e0a 100644 --- a/source/bin/restRoot.cxx +++ b/source/bin/restRoot.cxx @@ -25,6 +25,10 @@ int main(int argc, char* argv[]) { // set the env and debug status setenv("REST_VERSION", REST_RELEASE, 1); + printf("Setting verbose level to info. You may change level using `restRoot -v N`.\n"); + printf("Use `restRoot --help` for additional info.\n"); + gVerbose = StringToVerboseLevel("2"); + Bool_t loadMacros = false; for (int i = 1; i < argc; i++) { char* c = &argv[i][0]; @@ -100,6 +104,10 @@ int main(int argc, char* argv[]) { for (int i = 1; i < argc; i++) { const string opt = (string)argv[i]; + if (opt.at(0) == ('-') && opt.length() > 1 && opt.at(1) == ('-')) { + i++; + continue; + } if (opt.at(0) == ('-')) continue; if (opt.find("http") == string::npos && !TRestTools::fileExists(opt)) { diff --git a/source/framework/analysis/inc/TRestDataSetGainMap.h b/source/framework/analysis/inc/TRestDataSetGainMap.h new file mode 100644 index 000000000..50000a01b --- /dev/null +++ b/source/framework/analysis/inc/TRestDataSetGainMap.h @@ -0,0 +1,216 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) * + * For more information see https://gifna.unizar.es/trex * + * * + * REST is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * REST is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have a copy of the GNU General Public License along with * + * REST in $REST_PATH/LICENSE. * + * If not, see https://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +#ifndef REST_TRestDataSetGainMap +#define REST_TRestDataSetGainMap + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TRestDataSet.h" +#include "TRestMetadata.h" + +/// Metadata class to calculate,store and apply the gain corrected calibration of a group of detectors. +class TRestDataSetGainMap : public TRestMetadata { + public: + class Module; + + private: + std::string fObservable = ""; //"rawAna_ThresholdIntegral"; //< + std::string fSpatialObservableX = ""; //"hitsAna_xMean"; //< + std::string fSpatialObservableY = ""; //"hitsAna_yMean"; //< + + std::vector fModulesCal = {}; + std::string fCalibFileName = ""; + std::string fOutputFileName = ""; + + void Initialize() override; + void InitFromConfigFile() override; + + public: + std::set GetPlaneIDs() const; + std::set GetModuleIDs(const int planeId) const; + std::map> GetModuleIDs() const; + Int_t GetNumberOfPlanes() const { return GetPlaneIDs().size(); } + Int_t GetNumberOfModules() const { + int sum = 0; + for (auto pID : GetPlaneIDs()) sum += GetModuleIDs(pID).size(); + return sum; + } + + std::string GetCalibrationFileName() const { return fCalibFileName; } + std::string GetOutputFileName() const { return fOutputFileName; } + std::string GetObservable() const { return fObservable; } + std::string GetSpatialObservableX() const { return fSpatialObservableX; } + std::string GetSpatialObservableY() const { return fSpatialObservableY; } + + Module* GetModule(const int planeID, const int moduleID); + double GetSlopeParameter(const int planeID, const int moduleID, const double x, const double y); + double GetInterceptParameter(const int planeID, const int moduleID, const double x, const double y); + + void SetCalibrationFileName(const std::string& fileName) { fCalibFileName = fileName; } + void SetOutputFileName(const std::string& fileName) { fOutputFileName = fileName; } + void SetModuleCalibration(const Module& moduleCal); + void SetObservable(const std::string& observable) { fObservable = observable; } + void SetSpatialObservableX(const std::string& spatialObservableX) { + fSpatialObservableX = spatialObservableX; + } + void SetSpatialObservableY(const std::string& spatialObservableY) { + fSpatialObservableY = spatialObservableY; + } + + void Import(const std::string& fileName); + void Export(const std::string& fileName = ""); + + TRestDataSetGainMap& operator=(TRestDataSetGainMap& src); + + public: + void PrintMetadata() override; + + void GenerateGainMap(); + void CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName = ""); + + TRestDataSetGainMap(); + TRestDataSetGainMap(const char* configFilename, std::string name = ""); + ~TRestDataSetGainMap(); + + ClassDefOverride(TRestDataSetGainMap, 1); + + class Module { + private: + const TRestDataSetGainMap* p = nullptr; // fEnergyPeaks = {}; + std::vector fRangePeaks = {}; //{TVector2(230000, 650000), TVector2(40000, 230000)}; + TVector2 fCalibRange = TVector2(0, 0); //< // Calibration range + Int_t fNBins = 100; //< // Number of bins for the spectrum histograms + std::string fDefinitionCut = ""; //"TREXsides_tagId == 2"; //< + + Int_t fNumberOfSegmentsX = 1; //< + Int_t fNumberOfSegmentsY = 1; //< + TVector2 fReadoutRange = TVector2(-1, 246.24); //< // Readout dimensions + std::set fSplitX = {}; //< + std::set fSplitY = {}; //< + + std::string fDataSetFileName = ""; //< // File name for the calibration dataset + + std::vector> fSlope = {}; //< + std::vector> fIntercept = {}; //< + + bool fZeroPoint = false; //< Zero point will be automatically added if there are less than 2 peaks + bool fAutoRangePeaks = true; //< Automatic range peaks + std::vector> fSegSpectra = {}; + std::vector> fSegLinearFit = {}; + + public: + void AddPeak(const double& energyPeak, const TVector2& rangePeak = TVector2(0, 0)) { + fEnergyPeaks.push_back(energyPeak); + fRangePeaks.push_back(rangePeak); + } + + void LoadConfigFromTiXmlElement(const TiXmlElement* module); + + std::pair GetIndexMatrix(const double x, const double y) const; + double GetSlope(const double x, const double y) const; + double GetIntercept(const double x, const double y) const; + + Int_t GetPlaneId() const { return fPlaneId; } + Int_t GetModuleId() const { return fModuleId; } + std::string GetObservable() const { return p->fObservable; } + std::string GetSpatialObservableX() const { return p->fSpatialObservableX; } + std::string GetSpatialObservableY() const { return p->fSpatialObservableY; } + inline std::string GetModuleDefinitionCut() const { return fDefinitionCut; } + Int_t GetNumberOfSegmentsX() const { return fNumberOfSegmentsX; } + Int_t GetNumberOfSegmentsY() const { return fNumberOfSegmentsY; } + + std::set GetSplitX() const { return fSplitX; } + std::set GetSplitY() const { return fSplitY; } + std::string GetDataSetFileName() const { return fDataSetFileName; } + TVector2 GetReadoutRangeVar() const { return fReadoutRange; } + + void DrawSpectrum(const bool drawFits = true, const int color = -1, TCanvas* c = nullptr); + void DrawSpectrum(const TVector2& position, bool drawFits = true, int color = -1, + TCanvas* c = nullptr); + void DrawSpectrum(const int index_x, const int index_y, bool drawFits = true, int color = -1, + TCanvas* c = nullptr); + void DrawFullSpectrum(); + + void DrawLinearFit(); + void DrawLinearFit(const TVector2& position, TCanvas* c = nullptr); + void DrawLinearFit(const int index_x, const int index_y, TCanvas* c = nullptr); + + void DrawGainMap(const int peakNumber = 0); + + void Refit(const TVector2& position, const double energy, const TVector2& range); + void Refit(const size_t x, const size_t y, const size_t peakNumber, const TVector2& range); + void UpdateCalibrationFits(const size_t x, const size_t y); + + void SetPlaneId(const Int_t& planeId) { fPlaneId = planeId; } + void SetModuleId(const Int_t& moduleId) { fModuleId = moduleId; } + void SetModuleDefinitionCut(const std::string& moduleDefinitionCut) { + fDefinitionCut = moduleDefinitionCut; + } + void SetCalibrationRange(const TVector2& calibrationRange) { fCalibRange = calibrationRange; } + void SetNBins(const Int_t& nBins) { fNBins = nBins; } + void SetSplitX(); + void SetSplitY(); + void SetSplitX(const std::set& splitX); + void SetSplitY(const std::set& splitY); + void SetSplits(); + void SetSplits(const std::set& splitXandY) { + SetSplitX(splitXandY); + SetSplitY(splitXandY); + } + + void SetNumberOfSegmentsX(const Int_t& numberOfSegmentsX) { fNumberOfSegmentsX = numberOfSegmentsX; } + void SetNumberOfSegmentsY(const Int_t& numberOfSegmentsY) { fNumberOfSegmentsY = numberOfSegmentsY; } + + void SetDataSetFileName(const std::string& dataSetFileName) { fDataSetFileName = dataSetFileName; } + void SetReadoutRange(const TVector2& readoutRange) { fReadoutRange = readoutRange; } + void SetZeroPoint(const bool& ZeroPoint) { fZeroPoint = ZeroPoint; } + void SetAutoRangePeaks(const bool& autoRangePeaks) { fAutoRangePeaks = autoRangePeaks; } + + void Print() const; + + void GenerateGainMap(); + void Initialize(); + + Module() {} + Module(const TRestDataSetGainMap& parent) : p(&parent){}; + Module(const TRestDataSetGainMap& parent, const Int_t planeId, const Int_t moduleId) : p(&parent) { + SetPlaneId(planeId); + SetModuleId(moduleId); + }; + ~Module(){}; + }; +}; +#endif diff --git a/source/framework/analysis/inc/TRestDataSetOdds.h b/source/framework/analysis/inc/TRestDataSetOdds.h index 61210a3b9..7eed5fe1a 100644 --- a/source/framework/analysis/inc/TRestDataSetOdds.h +++ b/source/framework/analysis/inc/TRestDataSetOdds.h @@ -62,9 +62,18 @@ class TRestDataSetOdds : public TRestMetadata { void ComputeLogOdds(); + std::vector> GetOddsObservables(); + std::string GetOddsFile() { return fOddsFile; } + std::string GetDataSetName() { return fDataSetName; } + std::string GetOutputFileName() { return fOutputFileName; } + TRestCut* GetCut() { return fCut; } + inline void SetDataSetName(const std::string& dSName) { fDataSetName = dSName; } inline void SetOutputFileName(const std::string& outName) { fOutputFileName = outName; } inline void SetOddsFile(const std::string& oddsFile) { fOddsFile = oddsFile; } + inline void SetCut(TRestCut* cut) { fCut = cut; } + void SetOddsObservables(const std::vector>& obs); + void AddOddsObservable(const std::string& name, const TVector2& range, int nbins); TRestDataSetOdds(); TRestDataSetOdds(const char* configFilename, std::string name = ""); diff --git a/source/framework/analysis/src/TRestDataSetGainMap.cxx b/source/framework/analysis/src/TRestDataSetGainMap.cxx new file mode 100644 index 000000000..53a5fa0a0 --- /dev/null +++ b/source/framework/analysis/src/TRestDataSetGainMap.cxx @@ -0,0 +1,1285 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2016 GIFNA/TREX (University of Zaragoza) * + * For more information see https://gifna.unizar.es/trex * + * * + * REST is free software: you can redistribute it and/or modify * + * it under the terms of the GNU General Public License as published by * + * the Free Software Foundation, either version 3 of the License, or * + * (at your option) any later version. * + * * + * REST is distributed in the hope that it will be useful, * + * but WITHOUT ANY WARRANTY; without even the implied warranty of * + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * + * GNU General Public License for more details. * + * * + * You should have a copy of the GNU General Public License along with * + * REST in $REST_PATH/LICENSE. * + * If not, see https://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +///////////////////////////////////////////////////////////////////////// +/// TRestDataSetGainMap calculates and stores the calibration +/// parameters for a given detector with multiple (or just one) modules. +/// The modules are defined using the Module class (defined internally). +/// It performs a gain correction based on a spatial segmentation of the +/// detector module. This is useful forbig modules such as the ones used +/// in TREX-DM experiment. The input data files for this methods are +/// TRestDataSet for both calculating the calibration parameters and +/// performing the calibration of the desired events. +/// +/// To correct the gain inhomogeneities, the module is divided in +/// fNumberOfSegmentsX*fNumberOfSegmentsY segments. The energy peaks provided +/// are fitted in each segment. The events are associated to each segment based on +/// the observables fSpatialObservableX and fSpatialObservablesY. This results in +/// the following plot that can be obtain with the function +/// TRestDataSetGainMap::Module::DrawSpectrum() +/// +/// \htmlonly \endhtmlonly +/// ![Peak fitting of each segment. Plot obtain with +/// TRestDataSetGainMap::Module::DrawSpectrum()](drawSpectrum.png) +/// +/// Also, the peak position provides a gain map that can be plotted with the function +/// TRestDataSetGainMap::Module::DrawGainMap(peakNumber) +/// +/// \htmlonly \endhtmlonly +/// ![Gain map. Plot obtain with TRestDataSetGainMap::Module::DrawGainMap()](drawGainMap.png) +/// +/// The result is a better energy resolution with the gain corrected +/// calibration (red) than the plain calibration (blue). +/// +/// \htmlonly \endhtmlonly +/// ![Gain correction comparison.](gainCorrectionComparison.png) + +/// ### Parameters +/// * **calibFileName**: name of the file to use for the calibration. It should be a TRestDataSet +/// * **outputFileName**: name of the file to save this calibration metadata +/// * **observable**: name of the observable to be calibrated. It must be a branch of the calibration +/// TRestDataSet +/// * **spatialObservableX**: name of the observable to be used for the spatial segmentation along the X axis. +/// It must be a column (branch) of the calibration TRestDataSet +/// * **spatialObservableY**: name of the observable to be used for the spatial segmentation along the Y axis. +/// It must be a column (branch) of the calibration TRestDataSet +/// * **modulesCal**: vector of Module objects +/// +/// ### Examples +/// A example of RML definition of this class can be found inside the file +/// `REST_PATH/examples/calibrationCorrection.rml`. This have the following structure: +/// \code +/// +/// ... +/// +/// ... +/// +/// +/// ... +/// +/// +/// \endcode +/// +/// Example to calculate the calibration parameters over a calibration dataSet using restRoot: +/// \code +/// TRestDataSetGainMap gm ("calibrationCorrection.rml"); +/// gm.SetCalibrationFileName("myDataSet.root"); //if not already defined in rml file +/// gm.SetOutputFileName("myCalibration.root"); //if not already defined in rml file +/// gm.GenerateGainMap(); +/// gm.Export(); // gm.Export("anyOtherFileName.root") +/// \endcode +/// +/// Example to calibrate a dataSet with the previously calculated calibration parameters +/// and draw the result (using restRoot): +/// \code +/// TRestDataSetGainMap gm; +/// gm.Import("myCalibration.root"); +/// gm.CalibrateDataSet("dataSetToCalibrate.root", "calibratedDataSet.root"); +/// TRestDataSet ds; +/// ds.Import("calibratedDataSet.root"); +/// auto h = ds.GetDataFrame().Histo1D({"hname", "",100,-1,40.}, "calib_ThresholdIntegral"); +/// h->Draw(); +/// \endcode +/// +/// Example to refit manually the peaks of the gain map if any of them is not well fitted +/// (using restRoot): +/// \code +/// TRestDataSetGainMap gm; +/// gm.Import("myCalibration.root"); +/// // Draw all segmented spectra and check if any need a refit +/// for (auto pID : gm.GetPlaneIDs()) +/// for (auto mID : gm.GetModuleIDs(pID)) +/// gm.GetModule(pID,mID)->DrawSpectrum(); +/// // Draw only the desired spectrum (segment 0,0 in this case) +/// gm.GetModule(0,0)->DrawSpectrum(0, 0); +/// // Refit the desired peak (peak with energy 22.5 in this case) with a new range +/// TVector2 range(100000, 200000); // Define here the new range for the fit as you wish +/// gm.GetModule(0,0)->Refit(TVector2(0,0), 22.5, range); +/// // Check the result +/// gm.GetModule(0,0)->DrawSpectrum(TVector2(0.0,0.0)); // using x,y physical coord is possible +/// // Export the new calibration +/// gm.Export(); // gm.Export("anyOtherFileName.root") +/// \endcode +/// +///---------------------------------------------------------------------- +/// +/// REST-for-Physics - Software for Rare Event Searches Toolkit +/// +/// History of developments: +/// +/// 2023-September: First implementation of TRestDataSetGainMap +/// Álvaro Ezquerro +/// +/// \class TRestDataSetGainMap +/// \author: Álvaro Ezquerro aezquerro@unizar.es +/// +///
+/// + +#include "TRestDataSetGainMap.h" + +ClassImp(TRestDataSetGainMap); +/////////////////////////////////////////////// +/// \brief Default constructor +/// +TRestDataSetGainMap::TRestDataSetGainMap() { Initialize(); } + +///////////////////////////////////////////// +/// \brief Constructor loading data from a config file +/// +/// If no configuration path is defined using TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// The default behaviour is that the config file must be specified with +/// full path, absolute or relative. +/// +/// \param configFilename A const char* that defines the RML filename. +/// \param name The name of the metadata section. It will be used to find the +/// corresponding TRestDataSetGainMap section inside the RML. +/// +TRestDataSetGainMap::TRestDataSetGainMap(const char* configFilename, std::string name) + : TRestMetadata(configFilename) { + LoadConfigFromFile(fConfigFileName, name); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestDataSetGainMap::~TRestDataSetGainMap() {} + +void TRestDataSetGainMap::Initialize() { SetSectionName(this->ClassName()); } +/////////////////////////////////////////////// +/// \brief Initialization of TRestDataSetGainMap +/// members through a RML file +/// +void TRestDataSetGainMap::InitFromConfigFile() { + this->Initialize(); + TRestMetadata::InitFromConfigFile(); + + // Load Module from RML + TiXmlElement* moduleDefinition = GetElement("module"); + while (moduleDefinition != nullptr) { + fModulesCal.push_back(Module(*this)); + fModulesCal.back().LoadConfigFromTiXmlElement(moduleDefinition); + + moduleDefinition = GetNextElement(moduleDefinition); + } + + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata(); +} + +///////////////////////////////////////////// +/// \brief Function to calculate the calibration parameters of all modules +/// +void TRestDataSetGainMap::GenerateGainMap() { + for (auto& mod : fModulesCal) { + RESTInfo << "Calibrating plane " << mod.GetPlaneId() << " module " << mod.GetModuleId() << RESTendl; + mod.GenerateGainMap(); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) { + mod.DrawSpectrum(); + mod.DrawGainMap(); + } + } +} + +///////////////////////////////////////////// +/// \brief Function to calibrate a dataset. +/// +/// \param dataSetFileName the name of the root file where the TRestDataSet to be calibrated is stored. +/// \param outputFileName the name of the output (root) file where the calibrated TRestDataSet will be +/// exported. If empty, the output file will be named as the input file with the suffix "_cc". E.g. +/// "data/myDataSet.root" -> "data/myDataSet_cc.root". +/// +void TRestDataSetGainMap::CalibrateDataSet(const std::string& dataSetFileName, std::string outputFileName) { + if (fModulesCal.empty()) { + RESTError << "TRestDataSetGainMap::CalibrateDataSet: No modules defined." << RESTendl; + return; + } + + TRestDataSet dataSet; + dataSet.Import(dataSetFileName); + auto dataFrame = dataSet.GetDataFrame(); + + // Define a new column with the identifier (pmID) of the module for each row (event) + std::string pmIDname = (std::string)GetName() + "_pmID"; + std::string modCut = fModulesCal[0].GetModuleDefinitionCut(); + int pmID = fModulesCal[0].GetPlaneId() * 10 + fModulesCal[0].GetModuleId(); + + auto columnList = dataFrame.GetColumnNames(); + if (std::find(columnList.begin(), columnList.end(), pmIDname) == columnList.end()) + dataFrame = dataFrame.Define(pmIDname, modCut + " ? " + std::to_string(pmID) + " : -1"); + else + dataFrame = dataFrame.Redefine(pmIDname, modCut + " ? " + std::to_string(pmID) + " : -1"); + + for (size_t n = 1; n < fModulesCal.size(); n++) { + modCut = fModulesCal[n].GetModuleDefinitionCut(); + pmID = fModulesCal[n].GetPlaneId() * 10 + fModulesCal[n].GetModuleId(); + dataFrame = dataFrame.Redefine(pmIDname, (modCut + " ? " + std::to_string(pmID) + " : " + pmIDname)); + } + + // Define a new column with the calibrated observable + auto calibrate = [this](double val, double x, double y, int pmID) { + for (auto& m : fModulesCal) { + if (pmID == m.GetPlaneId() * 10 + m.GetModuleId()) + return m.GetSlope(x, y) * val + m.GetIntercept(x, y); + } + // RESTError << "TRestDataSetGainMap::CalibrateDataSet: Module not found for pmID " << pmID << + // RESTendl; + return std::numeric_limits::quiet_NaN(); + }; + std::string calibObsName = (std::string)GetName() + "_"; + calibObsName += GetObservable().erase(0, GetObservable().find("_") + 1); // remove the "rawAna_" part + dataFrame = dataFrame.Define(calibObsName, calibrate, + {fObservable, fSpatialObservableX, fSpatialObservableY, pmIDname}); + + dataSet.SetDataFrame(dataFrame); + + // Format the output file name and export the dataSet + if (outputFileName.empty()) outputFileName = dataSetFileName; + + if (outputFileName == dataSetFileName) { // TRestDataSet cannot be overwritten + outputFileName = outputFileName.substr(0, outputFileName.find_last_of(".")) + "_cc."; + outputFileName += TRestTools::GetFileNameExtension(dataSetFileName); + } + + RESTInfo << "Exporting calibrated dataSet to " << outputFileName << RESTendl; + dataSet.Export(outputFileName); + + // Add this TRestDataSetGainMap metadata to the output file + TFile* f = TFile::Open(outputFileName.c_str(), "UPDATE"); + this->Write(); + f->Close(); + delete f; +} + +///////////////////////////////////////////// +/// \brief Function to retrieve the module calibration with planeID and moduleID +/// +/// +TRestDataSetGainMap::Module* TRestDataSetGainMap::GetModule(const int planeID, const int moduleID) { + for (auto& i : fModulesCal) { + if (i.GetPlaneId() == planeID && i.GetModuleId() == moduleID) return &i; + } + RESTError << "No ModuleCalibration with planeID " << planeID << " and moduleID " << moduleID << RESTendl; + return nullptr; +} + +///////////////////////////////////////////// +/// \brief Function to get the slope parameter of the module with +/// planeID and moduleID at physical position (x,y) +/// +/// +double TRestDataSetGainMap::GetSlopeParameter(const int planeID, const int moduleID, const double x, + const double y) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetSlope(x, y); +} + +///////////////////////////////////////////// +/// \brief Function to get the intercept parameter of the module with +/// planeID and moduleID at physical position (x,y) +/// +/// +double TRestDataSetGainMap::GetInterceptParameter(const int planeID, const int moduleID, const double x, + const double y) { + Module* moduleCal = GetModule(planeID, moduleID); + if (moduleCal == nullptr) return 0; // return numeric_limits::quiet_NaN() + return moduleCal->GetIntercept(x, y); +} + +///////////////////////////////////////////// +/// \brief Function to get a list (set) of the plane IDs +/// +/// +std::set TRestDataSetGainMap::GetPlaneIDs() const { + std::set planeIDs; + for (const auto& mc : fModulesCal) planeIDs.insert(mc.GetPlaneId()); + return planeIDs; +} + +///////////////////////////////////////////// +/// \brief Function to get a list (set) of the module IDs +/// of the plane with planeID +/// +std::set TRestDataSetGainMap::GetModuleIDs(const int planeId) const { + std::set moduleIDs; + for (const auto& mc : fModulesCal) + if (mc.GetPlaneId() == planeId) moduleIDs.insert(mc.GetModuleId()); + return moduleIDs; +} + +///////////////////////////////////////////// +/// \brief Function to get the map of the module IDs for each plane ID +/// +std::map> TRestDataSetGainMap::GetModuleIDs() const { + std::map> moduleIds; + for (const int planeId : GetPlaneIDs()) + moduleIds.insert(std::pair>(planeId, GetModuleIDs(planeId))); + return moduleIds; +} + +TRestDataSetGainMap& TRestDataSetGainMap::operator=(TRestDataSetGainMap& src) { + SetName(src.GetName()); + fOutputFileName = src.GetOutputFileName(); + fObservable = src.GetObservable(); + fSpatialObservableX = src.GetSpatialObservableX(); + fSpatialObservableY = src.GetSpatialObservableY(); + fModulesCal.clear(); + for (auto pID : src.GetPlaneIDs()) + for (auto mID : src.GetModuleIDs(pID)) fModulesCal.push_back(*src.GetModule(pID, mID)); + return *this; +} + +///////////////////////////////////////////// +/// \brief Function to set a module calibration. If the module calibration +/// already exists (same planeId and moduleId), it will be replaced. +/// +void TRestDataSetGainMap::SetModuleCalibration(const Module& moduleCal) { + for (auto& i : fModulesCal) { + if (i.GetPlaneId() == moduleCal.GetPlaneId() && i.GetModuleId() == moduleCal.GetModuleId()) { + i = moduleCal; + return; + } + } + fModulesCal.push_back(moduleCal); +} + +///////////////////////////////////////////// +/// \brief Function to import the calibration parameters +/// from the root file fileName. +/// +void TRestDataSetGainMap::Import(const std::string& fileName) { + if (fileName.empty()) { + RESTError << "No input calibration file defined" << RESTendl; + return; + } + + if (TRestTools::isRootFile(fileName)) { + RESTInfo << "Opening " << fileName << RESTendl; + TFile* f = TFile::Open(fileName.c_str(), "READ"); + if (f == nullptr) { + RESTError << "Cannot open calibration file " << fileName << RESTendl; + return; + } + + TRestDataSetGainMap* cal = nullptr; + if (f != nullptr) { + TIter nextkey(f->GetListOfKeys()); + TKey* key; + while ((key = (TKey*)nextkey())) { + std::string kName = key->GetClassName(); + if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr && + REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSetGainMap")) { + cal = f->Get(key->GetName()); + *this = *cal; + } + } + } + } else + RESTError << "File extension not supported for " << fileName << RESTendl; + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Debug) PrintMetadata(); +} + +///////////////////////////////////////////// +/// \brief Function to export the calibration +/// to the file fileName. +/// \param fileName The name of the output file. If empty, the output file +/// will be the fOutputFileName class member. +/// +void TRestDataSetGainMap::Export(const std::string& fileName) { + if (!fileName.empty()) fOutputFileName = fileName; + if (fOutputFileName.empty()) { + RESTError << "No output file defined" << RESTendl; + return; + } + + if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") { + TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE"); + this->Write(GetName()); + f->Close(); + delete f; + RESTInfo << "Calibration saved to " << fOutputFileName << RESTendl; + } else + RESTError << "File extension for " << fOutputFileName << "is not supported." << RESTendl; +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the metadata members +/// +void TRestDataSetGainMap::PrintMetadata() { + TRestMetadata::PrintMetadata(); + RESTMetadata << " Calibration dataset: " << fCalibFileName << RESTendl; + RESTMetadata << " Output file: " << fOutputFileName << RESTendl; + RESTMetadata << " Number of planes: " << GetNumberOfPlanes() << RESTendl; + RESTMetadata << " Number of modules: " << GetNumberOfModules() << RESTendl; + RESTMetadata << " Calibration observable: " << fObservable << RESTendl; + RESTMetadata << " Spatial observable X: " << fSpatialObservableX << RESTendl; + RESTMetadata << " Spatial observable Y: " << fSpatialObservableY << RESTendl; + RESTMetadata << "-----------------------------------------------" << RESTendl; + for (auto& i : fModulesCal) i.Print(); + RESTMetadata << "***********************************************" << RESTendl; + RESTMetadata << RESTendl; +} + +///////////////////////////////////////////// +/// \brief Function to get the index of the matrix of calibration parameters +/// for a given x and y position on the detector plane. +/// +/// \param x A const double that defines the x position on the detector plane. +/// \param y A const double that defines the y position on the detector plane. +/// +std::pair TRestDataSetGainMap::Module::GetIndexMatrix(const double x, const double y) const { + int index_x = -1, index_y = -1; + + if (fSplitX.upper_bound(x) != fSplitX.end()) { + index_x = std::distance(fSplitX.begin(), fSplitX.upper_bound(x)) - 1; + if (index_x < 0) { + RESTWarning << "index_x < 0 for x = " << x << " and fSplitX[0]=" << *fSplitX.begin() + << p->RESTendl; + index_x = 0; + } + } else { + RESTWarning << "x is out of split for x = " << x << p->RESTendl; + index_x = fSplitX.size() - 2; + } + + if (fSplitY.upper_bound(y) != fSplitY.end()) { + index_y = std::distance(fSplitY.begin(), fSplitY.upper_bound(y)) - 1; + if (index_y < 0) { + RESTWarning << "index_y < 0 for y = " << y << " and fSplitY[0]=" << *fSplitY.begin() + << p->RESTendl; + index_y = 0; + } + } else { + RESTWarning << "y is out of split for y = " << y << p->RESTendl; + index_y = fSplitY.size() - 2; + } + + return std::make_pair(index_x, index_y); +} +///////////////////////////////////////////// +/// \brief Function to get the calibration parameter slope for a +/// given x and y position on the detector plane. +/// +/// \param x A const double that defines the x position on the detector plane. +/// \param y A const double that defines the y position on the detector plane. +/// +double TRestDataSetGainMap::Module::GetSlope(const double x, const double y) const { + auto [index_x, index_y] = GetIndexMatrix(x, y); + if (fSlope.empty()) { + RESTError << "Calibration slope matrix is empty. Returning 0" << p->RESTendl; + return 0; + } + + if (index_x > (int)fSlope.size() || index_y > (int)fSlope.at(0).size()) { + RESTError << "Index out of range. Returning 0" << p->RESTendl; + return 0; + } + + return fSlope[index_x][index_y]; +} + +///////////////////////////////////////////// +/// \brief Function to get the calibration parameter intercept for a +/// given x and y position on the detector plane. +/// +/// \param x A const double that defines the x position on the detector plane. +/// \param y A const double that defines the y position on the detector plane. +/// +double TRestDataSetGainMap::Module::GetIntercept(const double x, const double y) const { + auto [index_x, index_y] = GetIndexMatrix(x, y); + if (fIntercept.empty()) { + RESTError << "Calibration constant matrix is empty. Returning 0" << p->RESTendl; + return 0; + } + + if (index_x > (int)fIntercept.size() || index_y > (int)fIntercept.at(0).size()) { + RESTError << "Index out of range. Returning 0" << p->RESTendl; + return 0; + } + + return fIntercept[index_x][index_y]; +} + +///////////////////////////////////////////// +/// \brief Function to set the class members for segmentation of +/// the detector plane along the X and Y axis. +void TRestDataSetGainMap::Module::SetSplits() { + SetSplitX(); + SetSplitY(); +} +///////////////////////////////////////////// +/// \brief Function to set the class members for segmentation of +/// the detector plane along the X axis. +/// +/// It uses the number of segments and the readout range to define the +/// edges of the segments. +void TRestDataSetGainMap::Module::SetSplitX() { + if (fNumberOfSegmentsX < 1) { + RESTError << "SetSplitX: fNumberOfSegmentsX must be >=1." << p->RESTendl; + return; + } + std::set split; + for (int i = 0; i <= fNumberOfSegmentsX; i++) { // <= so that the last segment is included + double x = + fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsX) * i; + split.insert(x); + } + SetSplitX(split); +} + +void TRestDataSetGainMap::Module::SetSplitX(const std::set& splitX) { + if (splitX.size() < 2) { + RESTError << "SetSplitX: split size must be >=2 (start and end of range must be included)." + << p->RESTendl; + return; + } + if (!fSlope.empty()) + RESTWarning << "SetSplitX: changing split but current gain map and calibration paremeters correspond " + "to previous splitting. Use GenerateGainMap() to update them." + << p->RESTendl; + fSplitX = splitX; + fNumberOfSegmentsX = fSplitX.size() - 1; +} +///////////////////////////////////////////// +/// \brief Function to set the class members for segmentation of +/// the detector plane along the Y axis. +/// +/// It uses the number of segments and the readout range to define the +/// edges of the segments. +void TRestDataSetGainMap::Module::SetSplitY() { + if (fNumberOfSegmentsY < 1) { + RESTError << "SetSplitY: fNumberOfSegmentsY must be >=1." << p->RESTendl; + return; + } + std::set split; + for (int i = 0; i <= fNumberOfSegmentsY; i++) { // <= so that the last segment is included + double y = + fReadoutRange.X() + ((fReadoutRange.Y() - fReadoutRange.X()) / (float)fNumberOfSegmentsY) * i; + split.insert(y); + } + SetSplitY(split); +} + +void TRestDataSetGainMap::Module::SetSplitY(const std::set& splitY) { + if (splitY.size() < 2) { + RESTError << "SetSplitY: split size must be >=2 (start and end of range must be included)." + << p->RESTendl; + return; + } + if (!fSlope.empty()) + RESTWarning << "SetSplitY: changing split but current gain map and calibration paremeters correspond " + "to previous splitting. Use GenerateGainMap() to update them." + << p->RESTendl; + fSplitY = splitY; + fNumberOfSegmentsY = fSplitY.size() - 1; +} + +///////////////////////////////////////////// +/// \brief Function that calculates the calibration parameters for each segment +/// defined at fSplitX and fSplitY ang generate their spectra and gain map. +/// +/// It uses the data of the observable fObservable from the TRestDataSet +/// at fDataSetFileName (or fCalibFileName if first is empty). +/// The segmentation is given by the splits. +/// +/// Ranges for peak fitting follows this logic: +/// 1. If fRangePeaks is defined and fAutoRangePeaks is false: fRangePeaks is used. +/// 2. If fRangePeaks is defined and fAutoRangePeaks is true: the fitting range +/// is calculated by the peak position found by TSpectrum inside fRangePeaks. +/// 3. If fRangePeaks is not defined and fAutoRangePeaks is false: use the peak position found by TSpectrum +/// and the peak position of the next peak to define the range. +/// 4. If fRangePeaks is not defined and fAutoRangePeaks is true: same as 3. +/// +void TRestDataSetGainMap::Module::GenerateGainMap() { + //--- Initial checks and settings --- + std::string dsFileName = fDataSetFileName; + if (dsFileName.empty()) dsFileName = p->GetCalibrationFileName(); + if (dsFileName.empty()) { + RESTError << "No calibration file defined" << p->RESTendl; + return; + } + + if (!TRestTools::fileExists(dsFileName)) { + RESTError << "Calibration file " << dsFileName << " does not exist." << p->RESTendl; + return; + } + if (!TRestTools::isDataSet(dsFileName)) RESTWarning << dsFileName << " is not a dataset." << p->RESTendl; + TRestDataSet dataSet; + dataSet.Import(dsFileName); + fDataSetFileName = dsFileName; + + SetSplits(); + + if (fSplitX.empty()) SetSplitX(); + if (fSplitY.empty()) SetSplitY(); + + //--- Get the calibration range if not provided (default is 0,0) --- + if (fCalibRange.X() >= fCalibRange.Y()) { + // Get spectrum for this file + std::string cut = fDefinitionCut; + if (cut.empty()) cut = "1"; + auto histo = dataSet.GetDataFrame().Filter(cut).Histo1D({"temp", "", fNBins, 0, 0}, GetObservable()); + std::unique_ptr hpunt = std::unique_ptr(static_cast(histo->Clone())); + // double xMin = hpunt->GetXaxis()->GetXmin(); + double xMax = hpunt->GetXaxis()->GetXmax(); + + // Reduce the range to avoid the possible empty (nCounts<1%) end part of the spectrum + double fraction = 1, nAtEndSpc = 0, step = 0.66; + while (nAtEndSpc * 1. / hpunt->Integral() < 0.01 && fraction > 0.001) { + fraction *= step; + nAtEndSpc = hpunt->Integral(hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax() * fraction), + hpunt->FindFixBin(hpunt->GetXaxis()->GetXmax())); + } + xMax = hpunt->GetXaxis()->GetXmax() * fraction / + step; // previous step is the last one that nAtEndSpc<1% + // Set the calibration range if needed + // fCalibRange.SetX(xMin); + fCalibRange.SetY(xMax); + hpunt.reset(); // delete hpunt; + RESTDebug << "Calibration range (auto)set to (" << fCalibRange.X() << "," << fCalibRange.Y() << ")" + << p->RESTendl; + } + + //--- Definition of histogram matrix --- + std::vector> h(fNumberOfSegmentsX, std::vector(fNumberOfSegmentsY, nullptr)); + for (size_t i = 0; i < h.size(); i++) { + for (size_t j = 0; j < h.at(0).size(); j++) { + std::string name = "hSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + + std::to_string(i) + "_" + std::to_string(j); + h[i][j] = new TH1F(name.c_str(), "", fNBins, fCalibRange.X(), + fCalibRange.Y()); // h[column][row] equivalent to h[x][y] + } + } + std::vector> calParSlope(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); + std::vector> calParIntercept(fNumberOfSegmentsX, + std::vector(fNumberOfSegmentsY, 0)); + + // build the spectrum for each segment + auto itX = fSplitX.begin(); + for (size_t i = 0; i < h.size(); i++) { + auto itY = fSplitY.begin(); + for (size_t j = 0; j < h.at(0).size(); j++) { + // Get the segment limits from the splits + auto xLower = *itX; + auto xUpper = *std::next(itX); + auto yLower = *itY; + auto yUpper = *std::next(itY); + + std::string segment_cut = ""; + if (!GetSpatialObservableX().empty()) + segment_cut += GetSpatialObservableX() + ">=" + std::to_string(xLower) + "&&" + + GetSpatialObservableX() + "<" + std::to_string(xUpper); + if (!GetSpatialObservableY().empty()) + segment_cut += "&&" + GetSpatialObservableY() + ">=" + std::to_string(yLower) + "&&" + + GetSpatialObservableY() + "<" + std::to_string(yUpper); + if (!fDefinitionCut.empty()) segment_cut += "&&" + fDefinitionCut; + if (segment_cut.empty()) segment_cut = "1"; + RESTExtreme << "Segment[" << i << "][" << j << "] cut: " << segment_cut << p->RESTendl; + auto histo = dataSet.GetDataFrame() + .Filter(segment_cut) + .Histo1D({"temp", "", h[i][j]->GetNbinsX(), h[i][j]->GetXaxis()->GetXmin(), + h[i][j]->GetXaxis()->GetXmax()}, + GetObservable()); + std::unique_ptr hpunt = std::unique_ptr(static_cast(histo->Clone())); + h[i][j]->Add(hpunt.get()); + hpunt.reset(); // delete hpunt; + itY++; + } + itX++; + } + + //--- Fit every peak energy for every segment --- + fSegLinearFit = std::vector(h.size(), std::vector(h.at(0).size(), nullptr)); + for (size_t i = 0; i < h.size(); i++) { + for (size_t j = 0; j < h.at(0).size(); j++) { + RESTExtreme << "Segment[" << i << "][" << j << "]" << p->RESTendl; + // Search for peaks --> peakPos + std::unique_ptr s(new TSpectrum(2 * fEnergyPeaks.size() + 1)); + std::vector peakPos; + s->Search(h[i][j], 2, "goff", 0.1); + for (int k = 0; k < s->GetNPeaks(); k++) peakPos.push_back(s->GetPositionX()[k]); + std::sort(peakPos.begin(), peakPos.end(), std::greater()); + const double ratio = peakPos.size() == 0 + ? 1 + : peakPos.front() / fEnergyPeaks.front(); // to estimate peak position + + // Initialize graph for linear fit + std::shared_ptr gr; + gr = std::shared_ptr(new TGraph()); + gr->SetName("grFit"); + gr->SetTitle((";" + GetObservable() + ";energy").c_str()); + + // Fit every energy peak + int c = 0; + double mu = 0; + for (const auto& energy : fEnergyPeaks) { + RESTExtreme << "\t fitting energy " << DoubleToString(energy, "%g") << p->RESTendl; + // estimation of the peak position is between start and end + double pos = energy * ratio; + double start = pos * 0.8; + double end = pos * 1.2; + if (fRangePeaks.at(c).X() < fRangePeaks.at(c).Y()) { // if range is defined use it + start = fRangePeaks.at(c).X(); + end = fRangePeaks.at(c).Y(); + } + + do { + if (fAutoRangePeaks) { + if (peakPos.size() > 0) { + // Find the peak position that is between start and end + pos = peakPos.at(0); + while (!(start < pos && pos < end)) { + // if none of the peak position is + // between start and end, use the greater one. + if (pos == peakPos.back()) { + pos = peakPos.at(0); + break; + } + pos = *std::next(std::find(peakPos.begin(), peakPos.end(), + pos)); // get next peak position + } + peakPos.erase(std::find(peakPos.begin(), peakPos.end(), + pos)); // remove this peak position from the list + // final estimation of the peak range (idem fitting window) with this peak + // position pos + start = pos * 0.8; + end = pos * 1.2; + const double relDist = peakPos.size() > 0 ? (pos - peakPos.front()) / pos : 999; + if (relDist < 0.2) { // if the next peak is too close reduce the window width + start = pos * (1 - relDist / 2); + end = pos * (1 + relDist / 2); + } + } + } + + std::string name = "g" + std::to_string(c); + TF1* g = new TF1(name.c_str(), "gaus", start, end); + RESTExtreme << "\t\tat " << DoubleToString(pos, "%.3g") << ". Range(" + << DoubleToString(start, "%.3g") << ", " << DoubleToString(end, "%.3g") << ")" + << p->RESTendl; + + if (h[i][j]->GetFunction(name.c_str())) // remove previous fit + h[i][j]->GetListOfFunctions()->Remove(h[i][j]->GetFunction(name.c_str())); + + h[i][j]->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + mu = g->GetParameter(1); + RESTExtreme << "\t\tgaus mean " << DoubleToString(mu, "%g") << p->RESTendl; + } while (fAutoRangePeaks && peakPos.size() > 0 && + !(start < mu && mu < end)); // avoid small peaks on main peak tail + gr->SetPoint(c++, mu, energy); + } + s.reset(); // delete s; + + if (fZeroPoint) gr->SetPoint(c++, 0, 0); + while (gr->GetN() < 2) { // minimun 2 points needed for linear fit + gr->SetPoint(c++, 0, 0); + SetZeroPoint(true); + RESTDebug << "Not enough points for linear fit. Adding and setting zero point to true" + << p->RESTendl; + } + + // Linear fit + std::unique_ptr linearFit; + linearFit = std::unique_ptr(new TF1("linearFit", "pol1")); + gr->Fit("linearFit", "SQ"); // Q for quiet mode + + fSegLinearFit.at(i).at(j) = (TGraph*)gr->Clone(); + + const double slope = linearFit->GetParameter(1); + const double intercept = linearFit->GetParameter(0); + calParSlope.at(i).at(j) = slope; + calParIntercept.at(i).at(j) = intercept; + } + } + fSegSpectra = h; + fSlope = calParSlope; + fIntercept = calParIntercept; +} + +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for a given segment of the module. +/// +/// \param position position along X and Y axes at the detector module (in physical units). +/// \param energyPeak The energy of the peak to be fitted (in physical units). +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::Refit(const TVector2& position, const double energyPeak, + const TVector2& range) { + auto [index_x, index_y] = GetIndexMatrix(position.X(), position.Y()); + int peakNumber = -1; + for (size_t i = 0; i < fEnergyPeaks.size(); i++) + if (fEnergyPeaks.at(i) == energyPeak) { + peakNumber = i; + break; + } + if (peakNumber == -1) { + RESTError << "Energy " << energyPeak << " not found in the list of energy peaks" << p->RESTendl; + return; + } + Refit((size_t)index_x, (size_t)index_y, (size_t)peakNumber, range); +} + +///////////////////////////////////////////// +/// \brief Function to fit again manually a peak for a given segment of the module. The +/// calibration curve is updated after the fit. +/// +/// \param x index along X-axis of the corresponding segment. +/// \param y index along Y-axis of the corresponding segment. +/// \param peakNumber The index of the peak to be fitted. +/// \param range The range for the fitting of the peak (in the observables corresponding units). +/// +void TRestDataSetGainMap::Module::Refit(const size_t x, const size_t y, const size_t peakNumber, + const TVector2& range) { + if (fSegSpectra.empty()) { + RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; + return; + } + if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) { + RESTError << "Segment with index (" << x << ", " << y << ") not found" << p->RESTendl; + return; + } + if (peakNumber >= fEnergyPeaks.size()) { + RESTError << "Peak with index " << peakNumber << " not found" << p->RESTendl; + return; + } + + // Refit the desired peak + std::string name = "g" + std::to_string(peakNumber); + TF1* g = new TF1(name.c_str(), "gaus", range.X(), range.Y()); + TH1F* h = fSegSpectra.at(x).at(y); + while (h->GetFunction(name.c_str())) // clear previous fits for this peakNumber + h->GetListOfFunctions()->Remove(h->GetFunction(name.c_str())); + h->Fit(g, "R+Q0"); // use 0 to not draw the fit but save it + + // Change the point of the graph + UpdateCalibrationFits(x, y); +} + +///////////////////////////////////////////// +/// \brief Function to update the calibration curve for a given segment of the module. The calibration +/// curve is cleared and then the means of the gaussian fits for each energy peak are added. If there are +/// less than 2 fits, zero points are added. Then, the calibration curve is refitted (linearFit). +/// +/// \param x index along X-axis of the corresponding segment. +/// \param y index along Y-axis of the corresponding segment. +/// +void TRestDataSetGainMap::Module::UpdateCalibrationFits(const size_t x, const size_t y) { + if (fSegSpectra.empty()) { + RESTError << "No gain map found. Use GenerateGainMap() first." << p->RESTendl; + return; + } + if (x >= fSegSpectra.size() || y >= fSegSpectra.at(0).size()) { + RESTError << "Segment with index (" << x << ", " << y << ") not found" << p->RESTendl; + return; + } + + TH1F* h = fSegSpectra.at(x).at(y); + TGraph* gr = fSegLinearFit.at(x).at(y); + + // Clear the points of the graph + for (size_t i = 0; i < fEnergyPeaks.size(); i++) gr->RemovePoint(i); + // Add the new points to the graph + int c = 0; + for (size_t i = 0; i < fEnergyPeaks.size(); i++) { + std::string fitName = (std::string) "g" + std::to_string(i); + TF1* g = h->GetFunction(fitName.c_str()); + if (!g) { + RESTWarning << "No fit ( " << fitName << " ) found for energy peak " << fEnergyPeaks[i] + << " in segment " << x << "," << y << p->RESTendl; + continue; + } + gr->SetPoint(c++, g->GetParameter(1), fEnergyPeaks[i]); + } + + // Add zero points if needed (if there are less than 2 points) + while (gr->GetN() < 2) { + gr->SetPoint(c++, 0, 0); + RESTWarning << "Not enough points for linear fit at segment (" << x << ", " << y + << "). Adding zero point." << p->RESTendl; + } + + // Refit the calibration curve + TF1* lf = nullptr; + if (gr->GetFunction("linearFit")) + lf = gr->GetFunction("linearFit"); + else + lf = new TF1("linearFit", "pol1"); + gr->Fit(lf, "SQ"); // Q for quiet mode + fSlope.at(x).at(y) = lf->GetParameter(1); + fIntercept.at(x).at(y) = lf->GetParameter(0); +} + +///////////////////////////////////////////// +/// \brief Function to read the parameters from the RML element (TiXmlElement) +/// and set those class members. +/// +/// \param module A const TiXmlElement pointer that contains the information. +/// Example of this RML element: +/// +/// +/// +/// +void TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement(const TiXmlElement* module) { + if (module == nullptr) { + RESTError << "TRestDataSetGainMap::Module::LoadConfigFromTiXmlElement: module is nullptr" + << p->RESTendl; + return; + } + + std::string el = !module->Attribute("planeId") ? "Not defined" : module->Attribute("planeId"); + if (!(el.empty() || el == "Not defined")) this->SetPlaneId(StringToInteger(el)); + el = !module->Attribute("moduleId") ? "Not defined" : module->Attribute("moduleId"); + if (!(el.empty() || el == "Not defined")) this->SetModuleId(StringToInteger(el)); + + el = !module->Attribute("moduleDefinitionCut") ? "Not defined" : module->Attribute("moduleDefinitionCut"); + if (!(el.empty() || el == "Not defined")) this->SetModuleDefinitionCut(el); + + el = !module->Attribute("numberOfSegmentsX") ? "Not defined" : module->Attribute("numberOfSegmentsX"); + if (!(el.empty() || el == "Not defined")) this->SetNumberOfSegmentsX(StringToInteger(el)); + el = !module->Attribute("numberOfSegmentsY") ? "Not defined" : module->Attribute("numberOfSegmentsY"); + if (!(el.empty() || el == "Not defined")) this->SetNumberOfSegmentsY(StringToInteger(el)); + el = !module->Attribute("readoutRange") ? "Not defined" : module->Attribute("readoutRange"); + if (!(el.empty() || el == "Not defined")) this->SetReadoutRange(StringTo2DVector(el)); + + el = !module->Attribute("calibRange") ? "Not defined" : module->Attribute("calibRange"); + if (!(el.empty() || el == "Not defined")) this->SetCalibrationRange(StringTo2DVector(el)); + el = !module->Attribute("nBins") ? "Not defined" : module->Attribute("nBins"); + if (!(el.empty() || el == "Not defined")) this->SetNBins(StringToInteger(el)); + + el = !module->Attribute("dataSetFileName") ? "Not defined" : module->Attribute("dataSetFileName"); + if (!(el.empty() || el == "Not defined")) this->SetDataSetFileName(el); + + el = !module->Attribute("zeroPoint") ? "Not defined" : module->Attribute("zeroPoint"); + if (!(el.empty() || el == "Not defined")) this->SetZeroPoint(ToLower(el) == "true"); + el = !module->Attribute("autoRangePeaks") ? "Not defined" : module->Attribute("autoRangePeaks"); + if (!(el.empty() || el == "Not defined")) this->SetAutoRangePeaks(ToLower(el) == "true"); + + // Get peaks energy and range + TiXmlElement* peakDefinition = (TiXmlElement*)module->FirstChildElement("peak"); + while (peakDefinition != nullptr) { + double energy = 0; + TVector2 range = TVector2(0, 0); + + std::string ell = + !peakDefinition->Attribute("energy") ? "Not defined" : peakDefinition->Attribute("energy"); + if (ell.empty() || ell == "Not defined") { + RESTError << "< peak variable key does not contain energy!" << p->RESTendl; + exit(1); + } + energy = StringToDouble(ell); + + ell = !peakDefinition->Attribute("range") ? "Not defined" : peakDefinition->Attribute("range"); + if (!(ell.empty() || ell == "Not defined")) range = StringTo2DVector(ell); + + this->AddPeak(energy, range); + peakDefinition = (TiXmlElement*)peakDefinition->NextSiblingElement(); + } +} + +void TRestDataSetGainMap::Module::DrawSpectrum(const TVector2& position, bool drawFits, int color, + TCanvas* c) { + std::pair index = GetIndexMatrix(position.X(), position.Y()); + DrawSpectrum(index.first, index.second, drawFits, color, c); +} + +void TRestDataSetGainMap::Module::DrawSpectrum(const int index_x, const int index_y, bool drawFits, int color, + TCanvas* c) { + if (fSegSpectra.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + if (index_x < 0 || index_y < 0 || index_x >= (int)fSegSpectra.size() || + index_y >= (int)fSegSpectra.at(index_x).size()) { + RESTError << "Index out of range." << p->RESTendl; + return; + } + + if (!c) { + std::string t = "spectrum_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + + std::to_string(index_x) + "_" + std::to_string(index_y); + c = new TCanvas(t.c_str(), t.c_str()); + } + + auto xLower = *std::next(fSplitX.begin(), index_x); + auto xUpper = *std::next(fSplitX.begin(), index_x + 1); + auto yLower = *std::next(fSplitY.begin(), index_y); + auto yUpper = *std::next(fSplitY.begin(), index_y + 1); + std::string tH = "Spectrum x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") + + ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" + + GetObservable() + ";counts"; + fSegSpectra[index_x][index_y]->SetTitle(tH.c_str()); + + if (color > 0) fSegSpectra[index_x][index_y]->SetLineColor(color); + size_t colorT = fSegSpectra[index_x][index_y]->GetLineColor(); + fSegSpectra[index_x][index_y]->Draw("same"); + + if (drawFits) + for (size_t c = 0; c < fEnergyPeaks.size(); c++) { + auto fit = fSegSpectra[index_x][index_y]->GetFunction(("g" + std::to_string(c)).c_str()); + if (!fit) RESTWarning << "Fit for energy peak" << fEnergyPeaks[c] << " not found." << p->RESTendl; + if (!fit) continue; + fit->SetLineColor(c + 2 != colorT ? c + 2 : ++colorT); /* does not work with kRed, kBlue, etc. + as they are not defined with the same number as the first 10 basic colors. See + https://root.cern.ch/doc/master/classTColor.html#C01 and + https://root.cern.ch/doc/master/classTColor.html#C02 */ + fit->Draw("same"); + } +} + +///////////////////////////////////////////// +/// \brief Function to draw the spectrum for each segment of the module +/// on the same canvas. The canvas is divided in fNumberOfSegmentsX x fNumberOfSegmentsY +/// pads. The segments are drawn with the bottom left corner corresponding to the +/// minimun x and y values of the readout range and top right corner corresponding +/// to the maximum x and y values of the readout range. +/// Tip: define a canvas and use this same canvas along different calls to this function +/// to draw the spectra of different modules on the same canvas. +/// Example: +/// TCanvas *myCanvas = new TCanvas(); +/// module1->DrawSpectrum(false, kBlue, myCanvas); +/// module2->DrawSpectrum(false, kRed, myCanvas); +/// +/// \param drawFits A bool to also draw the fits or not. +/// \param color An int to set the color of the spectra. If negative, +/// the color of the spectra is not changed. +/// \param c A TCanvas pointer to draw the spectra. If none (nullptr) is given, +/// a new one is created. +/// +void TRestDataSetGainMap::Module::DrawSpectrum(const bool drawFits, const int color, TCanvas* c) { + if (fSegSpectra.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + if (!c) { + std::string t = "spectrum_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + c = new TCanvas(t.c_str(), t.c_str()); + } + + size_t nPads = 0; + for (const auto& object : *c->GetListOfPrimitives()) + if (object->InheritsFrom(TVirtualPad::Class())) ++nPads; + if (nPads != 0 && nPads != fSegSpectra.size() * fSegSpectra.at(0).size()) { + RESTError << "Canvas " << c->GetName() << " has " << nPads << " pads, but " + << fSegSpectra.size() * fSegSpectra.at(0).size() << " are needed." << p->RESTendl; + return; + } else if (nPads == 0) + c->Divide(fSegSpectra.size(), fSegSpectra.at(0).size()); + + for (size_t i = 0; i < fSegSpectra.size(); i++) { + for (size_t j = 0; j < fSegSpectra[i].size(); j++) { + c->cd(i + 1 + fSegSpectra[i].size() * j); + DrawSpectrum(i, fSegSpectra[i].size() - 1 - j, drawFits, color, c); + } + } +} +void TRestDataSetGainMap::Module::DrawFullSpectrum() { + if (fSegSpectra.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + TH1F* sumHist = + new TH1F("fullSpc", "", fSegSpectra[0][0]->GetNbinsX(), fSegSpectra[0][0]->GetXaxis()->GetXmin(), + fSegSpectra[0][0]->GetXaxis()->GetXmax()); + + sumHist->SetTitle(("Full spectrum;" + GetObservable() + ";counts").c_str()); + for (size_t i = 0; i < fSegSpectra.size(); i++) { + for (size_t j = 0; j < fSegSpectra.at(0).size(); j++) { + sumHist->Add(fSegSpectra[i][j]); + } + } + std::string t = "fullSpc_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + TCanvas* c = new TCanvas(t.c_str(), t.c_str()); + c->cd(); + sumHist->Draw(); +} + +void TRestDataSetGainMap::Module::DrawLinearFit(const TVector2& position, TCanvas* c) { + std::pair index = GetIndexMatrix(position.X(), position.Y()); + DrawLinearFit(index.first, index.second, c); +} + +void TRestDataSetGainMap::Module::DrawLinearFit(const int index_x, const int index_y, TCanvas* c) { + if (fSegLinearFit.size() == 0) { + RESTError << "Spectra matrix is empty." << p->RESTendl; + return; + } + if (index_x < 0 || index_y < 0 || index_x >= (int)fSegLinearFit.size() || + index_y >= (int)fSegLinearFit.at(index_x).size()) { + RESTError << "Index out of range." << p->RESTendl; + return; + } + if (!c) { + std::string t = "linearFit_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId) + "_" + + std::to_string(index_x) + "_" + std::to_string(index_y); + c = new TCanvas(t.c_str(), t.c_str()); + } + auto xLower = *std::next(fSplitX.begin(), index_x); + auto xUpper = *std::next(fSplitX.begin(), index_x + 1); + auto yLower = *std::next(fSplitY.begin(), index_y); + auto yUpper = *std::next(fSplitY.begin(), index_y + 1); + std::string tH = "Linear Fit x=[" + DoubleToString(xLower, "%g") + ", " + DoubleToString(xUpper, "%g") + + ") y=[" + DoubleToString(yLower, "%g") + ", " + DoubleToString(yUpper, "%g") + ");" + + GetObservable() + ";energy"; + fSegLinearFit[index_x][index_y]->SetTitle(tH.c_str()); + fSegLinearFit[index_x][index_y]->Draw("AL*"); +} + +void TRestDataSetGainMap::Module::DrawLinearFit() { + std::string t = "linearFits_" + std::to_string(fPlaneId) + "_" + std::to_string(fModuleId); + TCanvas* myCanvas = new TCanvas(t.c_str(), t.c_str()); + myCanvas->Divide(fSegLinearFit.size(), fSegLinearFit.at(0).size()); + for (size_t i = 0; i < fSegLinearFit.size(); i++) { + for (size_t j = 0; j < fSegLinearFit[i].size(); j++) { + myCanvas->cd(i + 1 + fSegLinearFit[i].size() * j); + // fSegLinearFit[i][j]->Draw("AL*"); + DrawLinearFit(i, fSegSpectra[i].size() - 1 - j, myCanvas); + } + } +} + +void TRestDataSetGainMap::Module::DrawGainMap(const int peakNumber) { + if (peakNumber < 0 || peakNumber >= (int)fEnergyPeaks.size()) { + RESTError << "Peak number out of range (peakNumber should be between 0 and " + << fEnergyPeaks.size() - 1 << " )" << p->RESTendl; + return; + } + double peakEnergy = fEnergyPeaks[peakNumber]; + std::string title = "Gain map for energy " + DoubleToString(peakEnergy, "%g") + ";" + + GetSpatialObservableX() + ";" + GetSpatialObservableY(); // + " keV"; + std::string t = "gainMap" + std::to_string(peakNumber) + "_" + std::to_string(fPlaneId) + "_" + + std::to_string(fModuleId); + TCanvas* gainMap = new TCanvas(t.c_str(), t.c_str()); + gainMap->cd(); + TH2F* hGainMap = new TH2F(("h" + t).c_str(), title.c_str(), fNumberOfSegmentsY, fReadoutRange.X(), + fReadoutRange.Y(), fNumberOfSegmentsX, fReadoutRange.X(), fReadoutRange.Y()); + + const double peakPosRef = + fSegLinearFit[(fNumberOfSegmentsX - 1) / 2][(fNumberOfSegmentsY - 1) / 2]->GetPointX(peakNumber); + + auto itX = fSplitX.begin(); + for (size_t i = 0; i < fSegLinearFit.size(); i++) { + auto itY = fSplitY.begin(); + for (size_t j = 0; j < fSegLinearFit.at(0).size(); j++) { + auto xLower = *itX; + auto xUpper = *std::next(itX); + auto yLower = *itY; + auto yUpper = *std::next(itY); + hGainMap->Fill((xUpper + xLower) / 2., (yUpper + yLower) / 2., + fSegLinearFit[i][j]->GetPointX(peakNumber) / peakPosRef); + itY++; + } + itX++; + } + hGainMap->SetStats(0); + hGainMap->Draw("colz"); + hGainMap->SetBarOffset(0.2); + hGainMap->Draw("TEXT SAME"); +} + +///////////////////////////////////////////// +/// \brief Prints on screen the information about the +/// members of Module +/// +void TRestDataSetGainMap::Module::Print() const { + RESTMetadata << "-----------------------------------------------" << p->RESTendl; + RESTMetadata << " Plane ID: " << fPlaneId << p->RESTendl; + RESTMetadata << " Module ID: " << fModuleId << p->RESTendl; + RESTMetadata << " Definition cut: " << fDefinitionCut << p->RESTendl; + RESTMetadata << p->RESTendl; + + RESTMetadata << " Calibration dataset: " << fDataSetFileName << p->RESTendl; + RESTMetadata << p->RESTendl; + + RESTMetadata << " Energy peaks: "; + for (const auto& peak : fEnergyPeaks) RESTMetadata << peak << ", "; + RESTMetadata << p->RESTendl; + bool anyRange = false; + for (const auto& r : fRangePeaks) + if (r.X() < r.Y()) { + RESTMetadata << " Range peaks: "; + anyRange = true; + break; + } + if (anyRange) + for (const auto& r : fRangePeaks) RESTMetadata << "(" << r.X() << ", " << r.Y() << ") "; + if (anyRange) RESTMetadata << p->RESTendl; + RESTMetadata << " Auto range peaks: " << (fAutoRangePeaks ? "true" : "false") << p->RESTendl; + RESTMetadata << " Zero point: " << (fZeroPoint ? "true" : "false") << p->RESTendl; + RESTMetadata << " Calibration range: (" << fCalibRange.X() << ", " << fCalibRange.Y() << " )" + << p->RESTendl; + RESTMetadata << " Number of bins: " << fNBins << p->RESTendl; + RESTMetadata << p->RESTendl; + + RESTMetadata << " Number of segments X: " << fNumberOfSegmentsX << p->RESTendl; + RESTMetadata << " Number of segments Y: " << fNumberOfSegmentsY << p->RESTendl; + // RESTMetadata << " Draw: " << (fDrawVar ? "true" : "false") << p->RESTendl; + RESTMetadata << " Readout range (" << fReadoutRange.X() << ", " << fReadoutRange.Y() << " )" + << p->RESTendl; + RESTMetadata << "SplitX: "; + for (auto& i : fSplitX) { + RESTMetadata << " " << i; + } + RESTMetadata << p->RESTendl; + RESTMetadata << "SplitY: "; + for (auto& i : fSplitY) { + RESTMetadata << " " << i; + } + RESTMetadata << p->RESTendl; + RESTMetadata << p->RESTendl; + + RESTMetadata << " Slope: " << p->RESTendl; + size_t maxSize = 0; + for (auto& x : fSlope) + if (maxSize < x.size()) maxSize = x.size(); + for (size_t j = 0; j < maxSize; j++) { + for (size_t k = 0; k < fSlope.size(); k++) { + if (j < fSlope[k].size()) + RESTMetadata << DoubleToString(fSlope[k][fSlope[k].size() - 1 - j], "%.3e") << " "; + else + RESTMetadata << " "; + } + RESTMetadata << p->RESTendl; + } + RESTMetadata << " Intercept: " << p->RESTendl; + maxSize = 0; + for (auto& x : fIntercept) + if (maxSize < x.size()) maxSize = x.size(); + for (size_t j = 0; j < maxSize; j++) { + for (size_t k = 0; k < fIntercept.size(); k++) { + if (j < fIntercept[k].size()) + RESTMetadata << DoubleToString(fIntercept[k][fIntercept[k].size() - 1 - j], "%+.3e") << " "; + else + RESTMetadata << " "; + } + RESTMetadata << p->RESTendl; + } + RESTMetadata << "-----------------------------------------------" << p->RESTendl; +} diff --git a/source/framework/analysis/src/TRestDataSetOdds.cxx b/source/framework/analysis/src/TRestDataSetOdds.cxx index e51349dcc..f0d59bb99 100644 --- a/source/framework/analysis/src/TRestDataSetOdds.cxx +++ b/source/framework/analysis/src/TRestDataSetOdds.cxx @@ -209,14 +209,18 @@ void TRestDataSetOdds::ComputeLogOdds() { if (fOddsFile.empty()) { auto DF = dataSet.MakeCut(fCut); + RESTInfo << "Generating PDFs for dataset: " << fDataSetName << RESTendl; for (size_t i = 0; i < fObsName.size(); i++) { const std::string obsName = fObsName[i]; const TVector2 range = fObsRange[i]; const std::string histName = "h" + obsName; const int nBins = fObsNbins[i]; + RESTDebug << "\tGenerating PDF for " << obsName << " with range: (" << range.X() << ", " + << range.Y() << ") and nBins: " << nBins << RESTendl; auto histo = DF.Histo1D({histName.c_str(), histName.c_str(), nBins, range.X(), range.Y()}, obsName); TH1F* h = static_cast(histo->DrawClone()); + RESTDebug << "\tNormalizing by integral = " << h->Integral() << RESTendl; h->Scale(1. / h->Integral()); fHistos[obsName] = h; } @@ -226,7 +230,7 @@ void TRestDataSetOdds::ComputeLogOdds() { RESTError << "Cannot open calibration odds file " << fOddsFile << RESTendl; exit(1); } - std::cout << "Opening " << fOddsFile << std::endl; + RESTInfo << "Opening " << fOddsFile << " as oddsFile." << RESTendl; for (size_t i = 0; i < fObsName.size(); i++) { const std::string obsName = fObsName[i]; const std::string histName = "h" + obsName; @@ -237,6 +241,7 @@ void TRestDataSetOdds::ComputeLogOdds() { auto df = dataSet.GetDataFrame(); std::string totName = ""; + RESTDebug << "Computing log odds from " << fDataSetName << RESTendl; for (const auto& [obsName, histo] : fHistos) { const std::string oddsName = "odds_" + obsName; auto GetLogOdds = [&histo = histo](double val) { @@ -244,6 +249,12 @@ void TRestDataSetOdds::ComputeLogOdds() { if (odds == 0) return 1000.; return log(1. - odds) - log(odds); }; + + if (df.GetColumnType(obsName) != "Double_t") { + RESTWarning << "Column " << obsName << " is not of type 'double'. It will be converted." + << RESTendl; + df = df.Redefine(obsName, "static_cast(" + obsName + ")"); + } df = df.Define(oddsName, GetLogOdds, {obsName}); auto h = df.Histo1D(oddsName); @@ -251,27 +262,60 @@ void TRestDataSetOdds::ComputeLogOdds() { totName += oddsName; } + RESTDebug << "Computing total log odds" << RESTendl; + RESTDebug << "\tTotal log odds = " << totName << RESTendl; df = df.Define("odds_total", totName); dataSet.SetDataFrame(df); if (!fOutputFileName.empty()) { if (TRestTools::GetFileNameExtension(fOutputFileName) == "root") { + RESTDebug << "Exporting dataset to " << fOutputFileName << RESTendl; dataSet.Export(fOutputFileName); TFile* f = TFile::Open(fOutputFileName.c_str(), "UPDATE"); this->Write(); + RESTDebug << "Writing histograms to " << fOutputFileName << RESTendl; for (const auto& [obsName, histo] : fHistos) histo->Write(); f->Close(); } } } +std::vector> TRestDataSetOdds::GetOddsObservables() { + std::vector> obs; + for (size_t i = 0; i < fObsName.size(); i++) { + if (i >= fObsName.size() || i >= fObsRange.size() || i >= fObsNbins.size()) { + RESTError << "Sizes for observables names, ranges and bins do not match!" << RESTendl; + break; + } + obs.push_back(std::make_tuple(fObsName[i], fObsRange[i], fObsNbins[i])); + } + return obs; +} + +void TRestDataSetOdds::AddOddsObservable(const std::string& name, const TVector2& range, int nbins) { + fObsName.push_back(name); + fObsRange.push_back(range); + fObsNbins.push_back(nbins); +} + +void TRestDataSetOdds::SetOddsObservables(const std::vector>& obs) { + fObsName.clear(); + fObsRange.clear(); + fObsNbins.clear(); + for (const auto& [name, range, nbins] : obs) AddOddsObservable(name, range, nbins); +} + ///////////////////////////////////////////// /// \brief Prints on screen the information about the metadata members of TRestDataSetOdds /// void TRestDataSetOdds::PrintMetadata() { TRestMetadata::PrintMetadata(); + // if (fCut) fCut->PrintMetadata(); + if (!fOddsFile.empty()) RESTMetadata << " Odds file: " << fOddsFile << RESTendl; + RESTMetadata << " DataSet file: " << fDataSetName << RESTendl; + RESTMetadata << " Observables to compute: " << RESTendl; for (size_t i = 0; i < fObsName.size(); i++) { RESTMetadata << fObsName[i] << "; Range: (" << fObsRange[i].X() << ", " << fObsRange[i].Y() diff --git a/source/framework/core/inc/TRestDataSet.h b/source/framework/core/inc/TRestDataSet.h index 94f0a490d..4efe5699e 100644 --- a/source/framework/core/inc/TRestDataSet.h +++ b/source/framework/core/inc/TRestDataSet.h @@ -105,6 +105,10 @@ class TRestDataSet : public TRestMetadata { /// A list of new columns together with its corresponding expressions added to the dataset std::vector> fColumnNameExpressions; + /// A flag to enable Multithreading during dataframe generation + Bool_t fMT = false; //< + + inline auto GetAddedColumns() const { return fColumnNameExpressions; } /// The resulting RDF::RNode object after initialization ROOT::RDF::RNode fDataSet = ROOT::RDataFrame(0); //! @@ -125,6 +129,8 @@ class TRestDataSet : public TRestMetadata { void SetDataFrame(const ROOT::RDF::RNode& dS) { fDataSet = dS; } + void EnableMultiThreading(Bool_t enable = true) { fMT = enable; } + /// Gives access to the tree TTree* GetTree() const { if (fTree == nullptr) { @@ -169,6 +175,7 @@ class TRestDataSet : public TRestMetadata { inline void SetQuantity(const std::map& quantity) { fQuantity = quantity; } TRestDataSet& operator=(TRestDataSet& dS); + Bool_t Merge(const TRestDataSet& dS); void Import(const std::string& fileName); void Import(std::vector fileNames); void Export(const std::string& filename); diff --git a/source/framework/core/inc/TRestMetadata.h b/source/framework/core/inc/TRestMetadata.h index e1e1910ef..d87a60139 100644 --- a/source/framework/core/inc/TRestMetadata.h +++ b/source/framework/core/inc/TRestMetadata.h @@ -72,7 +72,7 @@ class TRestMetadata : public TNamed { void ReadEnvInElement(TiXmlElement* e, bool overwrite = true); void ReadElement(TiXmlElement* e, bool recursive = false); void ReplaceForLoopVars(TiXmlElement* e, std::map forLoopVar); - void ExpandForLoopOnce(TiXmlElement* e, std::map forLoopVar); + void ExpandForLoopOnce(TiXmlElement* e, const std::map& forLoopVar); void ExpandForLoops(TiXmlElement* e, std::map forLoopVar); void ExpandIfSections(TiXmlElement* e); void ExpandIncludeFile(TiXmlElement* e); @@ -216,10 +216,10 @@ class TRestMetadata : public TNamed { void AddLog(std::string log = "", bool print = true); /// A metadata class may use this method to signal that something went wrong - void SetError(std::string message = "", bool print = true); + void SetError(std::string message = "", bool print = true, int maxPrint = 5); /// A metadata class may use this method to signal that something went wrong - void SetWarning(std::string message = "", bool print = true); + void SetWarning(std::string message = "", bool print = true, int maxPrint = 5); /// Returns a std::string containing the error message TString GetErrorMessage(); @@ -239,6 +239,10 @@ class TRestMetadata : public TNamed { TRestMetadata* InstantiateChildMetadata(int index, std::string pattern = ""); TRestMetadata* InstantiateChildMetadata(std::string pattern = "", std::string name = ""); + /// Merge the metadata information from another metadata object. + /// Needs to be implemented in the derived class. + virtual void Merge(const TRestMetadata&); + /// Making default settings. virtual void Initialize() {} diff --git a/source/framework/core/inc/TRestProcessRunner.h b/source/framework/core/inc/TRestProcessRunner.h index 66e09fef9..bb8ef5ab6 100644 --- a/source/framework/core/inc/TRestProcessRunner.h +++ b/source/framework/core/inc/TRestProcessRunner.h @@ -87,7 +87,7 @@ class TRestProcessRunner : public TRestMetadata { EndOfInit(); } void BeginOfInit(); - Int_t ReadConfig(std::string keydeclare, TiXmlElement* e); + Int_t ReadConfig(const std::string& keydeclare, TiXmlElement* e); void EndOfInit(); void PrintMetadata() override; diff --git a/source/framework/core/src/TRestDataSet.cxx b/source/framework/core/src/TRestDataSet.cxx index 4e7d0aa56..a3d13b94d 100644 --- a/source/framework/core/src/TRestDataSet.cxx +++ b/source/framework/core/src/TRestDataSet.cxx @@ -283,8 +283,6 @@ TRestDataSet::TRestDataSet() { Initialize(); } /// TRestDataSet::TRestDataSet(const char* cfgFileName, const std::string& name) : TRestMetadata(cfgFileName) { LoadConfigFromFile(fConfigFileName, name); - - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) PrintMetadata(); } /////////////////////////////////////////////// @@ -343,7 +341,10 @@ void TRestDataSet::GenerateDataSet() { std::sort(finalList.begin(), finalList.end()); finalList.erase(std::unique(finalList.begin(), finalList.end()), finalList.end()); - ROOT::EnableImplicitMT(); + if (fMT) + ROOT::EnableImplicitMT(); + else + ROOT::DisableImplicitMT(); RESTInfo << "Initializing dataset" << RESTendl; fDataSet = ROOT::RDataFrame("AnalysisTree", fFileSelection); @@ -621,8 +622,11 @@ void TRestDataSet::PrintMetadata() { if (!fColumnNameExpressions.empty()) { RESTMetadata << " New columns added to generated dataframe: " << RESTendl; RESTMetadata << " ---------------------------------------- " << RESTendl; - for (const auto& [cName, cExpression] : fColumnNameExpressions) - RESTMetadata << " - Name : " << cName << " Expression: " << cExpression << RESTendl; + for (const auto& [cName, cExpression] : fColumnNameExpressions) { + RESTMetadata << " - Name : " << cName << RESTendl; + RESTMetadata << " - Expression: " << cExpression << RESTendl; + RESTMetadata << " " << RESTendl; + } } if (fMergedDataset) { @@ -638,6 +642,12 @@ void TRestDataSet::PrintMetadata() { for (const auto& fn : fImportedFiles) RESTMetadata << " - " << fn << RESTendl; } + RESTMetadata << " " << RESTendl; + if (fMT) + RESTMetadata << " - Multithreading was enabled" << RESTendl; + else + RESTMetadata << " - Multithreading was NOT enabled" << RESTendl; + RESTMetadata << "----" << RESTendl; } @@ -894,12 +904,35 @@ TRestDataSet& TRestDataSet::operator=(TRestDataSet& dS) { fFilterEqualsTo = dS.GetFilterEqualsTo(); fQuantity = dS.GetQuantity(); fTotalDuration = dS.GetTotalTimeInSeconds(); - fFileSelection = dS.GetFileSelection(); fCut = dS.GetCut(); return *this; } +/////////////////////////////////////////////// +/// \brief This function merge different TRestDataSet +/// metadata in current dataSet +/// +Bool_t TRestDataSet::Merge(const TRestDataSet& dS) { + auto obsNames = GetObservablesList(); + for (const auto& obs : fObservablesList) { + if (std::find(obsNames.begin(), obsNames.end(), obs) != obsNames.end()) { + RESTError << "Cannot merge dataSets with different observable list " << RESTendl; + return false; + } + } + + if (fStartTime > dS.GetStartTime()) fStartTime = dS.GetStartTime(); + if (fEndTime < dS.GetEndTime()) fEndTime = dS.GetEndTime(); + + auto fileSelection = dS.GetFileSelection(); + fFileSelection.insert(fFileSelection.end(), fileSelection.begin(), fileSelection.end()); + + fTotalDuration += dS.GetTotalTimeInSeconds(); + + return true; +} + /////////////////////////////////////////////// /// \brief This function imports metadata from a root file /// it import metadata info from the previous dataSet @@ -921,8 +954,6 @@ void TRestDataSet::Import(const std::string& fileName) { if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr && REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSet")) { dS = file->Get(key->GetName()); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) - dS->PrintMetadata(); *this = *dS; } } @@ -933,7 +964,10 @@ void TRestDataSet::Import(const std::string& fileName) { return; } - ROOT::EnableImplicitMT(); + if (fMT) + ROOT::EnableImplicitMT(); + else + ROOT::DisableImplicitMT(); fDataSet = ROOT::RDataFrame("AnalysisTree", fileName); @@ -956,25 +990,48 @@ void TRestDataSet::Import(std::vector fileNames) { return; } - if (fileNames.size() == 0) return; - - TFile* file = TFile::Open(fileNames[0].c_str(), "READ"); - if (file != nullptr) { - TIter nextkey(file->GetListOfKeys()); - TKey* key; - while ((key = (TKey*)nextkey())) { - std::string kName = key->GetClassName(); - if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr && - REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSet")) { - TRestDataSet* dS = file->Get(key->GetName()); - if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) - dS->PrintMetadata(); - *this = *dS; + int count = 0; + auto it = fileNames.begin(); + while (it != fileNames.end()) { + std::string fileName = *it; + TFile* file = TFile::Open(fileName.c_str(), "READ"); + bool isValid = false; + if (file != nullptr) { + TIter nextkey(file->GetListOfKeys()); + TKey* key; + while ((key = (TKey*)nextkey())) { + std::string kName = key->GetClassName(); + if (REST_Reflection::GetClassQuick(kName.c_str()) != nullptr && + REST_Reflection::GetClassQuick(kName.c_str())->InheritsFrom("TRestDataSet")) { + TRestDataSet* dS = file->Get(key->GetName()); + if (GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Info) + dS->PrintMetadata(); + + if (count == 0) { + *this = *dS; + isValid = true; + } else { + isValid = Merge(*dS); + } + + if (isValid) count++; + } } + } else { + RESTError << "Cannot open " << fileName << RESTendl; } - } else { - RESTError << "Cannot open " << fileNames[0] << RESTendl; - exit(1); + + if (!isValid) { + RESTError << fileName << " is not a valid dataSet skipping..." << RESTendl; + it = fileNames.erase(it); + } else { + ++it; + } + } + + if (fileNames.empty()) { + RESTError << "File selection is empty, dataSet will not be imported " << RESTendl; + return; } RESTInfo << "Opening list of files. First file: " << fileNames[0] << RESTendl; diff --git a/source/framework/core/src/TRestMetadata.cxx b/source/framework/core/src/TRestMetadata.cxx index 4fcccc1fe..718f62daa 100644 --- a/source/framework/core/src/TRestMetadata.cxx +++ b/source/framework/core/src/TRestMetadata.cxx @@ -1065,8 +1065,10 @@ void TRestMetadata::ExpandIfSections(TiXmlElement* e) { /////////////////////////////////////////////// /// \brief Helper method for TRestMetadata::ExpandForLoops(). -void TRestMetadata::ExpandForLoopOnce(TiXmlElement* e, map forLoopVar) { - if (e == nullptr) return; +void TRestMetadata::ExpandForLoopOnce(TiXmlElement* e, const map& forLoopVar) { + if (e == nullptr) { + return; + } TiXmlElement* parele = (TiXmlElement*)e->Parent(); TiXmlElement* contentelement = e->FirstChildElement(); @@ -1098,7 +1100,7 @@ void TRestMetadata::ReplaceForLoopVars(TiXmlElement* e, map forL if (e == nullptr) return; RESTDebug << "Entering ... TRestMetadata::ReplaceForLoopVars" << RESTendl; - std::string parName = ""; + std::string parName; TiXmlAttribute* attr = e->FirstAttribute(); while (attr != nullptr) { const char* val = attr->Value(); @@ -1111,7 +1113,7 @@ void TRestMetadata::ReplaceForLoopVars(TiXmlElement* e, map forL if (strcmp(name, "name") != 0) { string outputBuffer = val; - if (outputBuffer.find("[") != string::npos || outputBuffer.find("]") != string::npos) { + if (outputBuffer.find('[') != string::npos || outputBuffer.find(']') != string::npos) { RESTError << "TRestMetadata::ReplaceForLoopVars. Old for-loop construction identified" << RESTendl; RESTError << "Please, replace [] variable nomenclature by ${}." << RESTendl; @@ -1134,7 +1136,7 @@ void TRestMetadata::ReplaceForLoopVars(TiXmlElement* e, map forL string proenv = forLoopVar.count(expression) > 0 ? forLoopVar[expression] : ""; - if (proenv != "") { + if (!proenv.empty()) { outputBuffer.replace(replacePos, replaceLen, proenv); endPosition = 0; } else { @@ -1203,11 +1205,11 @@ void TRestMetadata::ExpandForLoops(TiXmlElement* e, map forloopv if (fVerboseLevel >= TRestStringOutput::REST_Verbose_Level::REST_Extreme) parele->Print(stdout, 0); RESTDebug << "----end of for loop----" << RESTendl; - } else if (_in != "") { + } else if (!_in.empty()) { vector loopvars = Split(_in, ":"); RESTDebug << "----expanding for loop----" << RESTendl; - for (string loopvar : loopvars) { + for (const string& loopvar : loopvars) { forloopvar[_name] = loopvar; fVariables[_name] = loopvar; ExpandForLoopOnce(e, forloopvar); @@ -1250,18 +1252,22 @@ void TRestMetadata::ExpandIncludeFile(TiXmlElement* e) { if (_filename == nullptr) return; string filename; - if (string(_filename) == "server") { + if (string(_filename) == "server" || TRestTools::isURL(_filename)) { // Let TRestRun to retrieve data according to run number later-on - // if ((string) this->ClassName() == "TRestRun") return; // match the database, runNumber=0(default data), type="META_RML", tag=
auto url = gDataBase->query_data(DBEntry(0, "META_RML", e->Value())).value; + if (url.empty()) { + // don't really understand this "database" code, this just works + url = _filename; + } + filename = TRestTools::DownloadRemoteFile(url); } else { filename = SearchFile(_filename); } - if (filename == "") { + if (filename.empty()) { RESTError << "TRestMetadata::ExpandIncludeFile. Include file \"" << _filename << "\" does not exist!" << RESTendl; exit(1); @@ -1369,7 +1375,7 @@ void TRestMetadata::ExpandIncludeFile(TiXmlElement* e) { } ele = ele->NextSiblingElement(); } - // more than 1 elements found + // more than 1 element found if (eles.size() > 1) { RESTWarning << "(expand include file): find multiple xml sections with same name!" << RESTendl; @@ -2626,23 +2632,23 @@ void TRestMetadata::AddLog(string log, bool print) { } } -void TRestMetadata::SetError(string message, bool print) { +void TRestMetadata::SetError(string message, bool print, int maxPrint) { fError = true; fNErrors++; if (message != "") { fErrorMessage += message + "\n"; - if (print) { + if (print && fNErrors < maxPrint) { cout << message << endl; } } } -void TRestMetadata::SetWarning(string message, bool print) { +void TRestMetadata::SetWarning(string message, bool print, int maxPrint) { fWarning = true; fNWarnings++; if (message != "") { fWarningMessage += message + "\n"; - if (print) { + if (print && fNWarnings < maxPrint) { RESTWarning << message << RESTendl; } } @@ -2661,3 +2667,14 @@ TString TRestMetadata::GetWarningMessage() { else return "No warning!"; } + +void TRestMetadata::Merge(const TRestMetadata& metadata) { + if (!metadata.InheritsFrom(ClassName())) { + RESTError << "TRestMetadata::Merge. Metadata is not of type " << ClassName() << RESTendl; + exit(1); + } + + if (fName.IsNull()) { + fName = metadata.GetName(); + } +} diff --git a/source/framework/core/src/TRestProcessRunner.cxx b/source/framework/core/src/TRestProcessRunner.cxx index 55b7f1e31..85a7efbd0 100644 --- a/source/framework/core/src/TRestProcessRunner.cxx +++ b/source/framework/core/src/TRestProcessRunner.cxx @@ -22,6 +22,8 @@ ///
////////////////////////////////////////////////////////////////////////// +#include + #include "Math/MinimizerOptions.h" #include "TBranchRef.h" #include "TInterpreter.h" @@ -169,7 +171,7 @@ void TRestProcessRunner::BeginOfInit() { // fUsePauseMenu = StringToBool(GetParameter("usePauseMenu", "OFF")); if (!fUsePauseMenu || fVerboseLevel >= TRestStringOutput::REST_Verbose_Level::REST_Debug) fProcStatus = kIgnore; - if (fOutputAnalysisStorage == false) { + if (!fOutputAnalysisStorage) { RESTError << "output analysis must be turned on to process data!" << RESTendl; exit(1); } @@ -183,8 +185,12 @@ void TRestProcessRunner::BeginOfInit() { // fOutputItem = Split(GetParameter("treeBranches", // "inputevent:outputevent:inputanalysis"), ":"); - if (fThreadNumber < 1) fThreadNumber = 1; - if (fThreadNumber > 15) fThreadNumber = 15; + if (fThreadNumber < 1) { + fThreadNumber = 1; + } + if (fThreadNumber > 15) { + fThreadNumber = 15; + } for (int i = 0; i < fThreadNumber; i++) { TRestThread* t = new TRestThread(); @@ -204,7 +210,7 @@ void TRestProcessRunner::BeginOfInit() { /// InstantiateProcess() The processes will be added into each TRestThread /// instance. If the process is external process, then it will be sent to /// TRestRun. -Int_t TRestProcessRunner::ReadConfig(string keydeclare, TiXmlElement* e) { +Int_t TRestProcessRunner::ReadConfig(const string& keydeclare, TiXmlElement* e) { if (keydeclare == "addProcess") { string active = GetParameter("value", e, ""); if (active != "" && ToUpper(active) != "ON") return 0; @@ -388,9 +394,9 @@ void TRestProcessRunner::RunProcess() { fEventTree->SetName("EventTree"); fEventTree->SetTitle("REST Event Tree"); fEventTree->SetDirectory(fOutputDataFile); - fEventTree->SetMaxTreeSize(100000000000LL > fFileSplitSize * 2 - ? 100000000000LL - : fFileSplitSize * 2); // the default size is 100GB + TTree::SetMaxTreeSize(100000000000LL > fFileSplitSize * 2 + ? 100000000000LL + : fFileSplitSize * 2); // the default size is 100GB } else { fEventTree = nullptr; } @@ -398,7 +404,8 @@ void TRestProcessRunner::RunProcess() { // initialize analysis tree fAnalysisTree = new TRestAnalysisTree("AnalysisTree", "REST Process Analysis Tree"); fAnalysisTree->SetDirectory(fOutputDataFile); - fAnalysisTree->SetMaxTreeSize(100000000000LL > fFileSplitSize * 2 ? 100000000000LL : fFileSplitSize * 2); + TRestAnalysisTree::SetMaxTreeSize(100000000000LL > fFileSplitSize * 2 ? 100000000000LL + : fFileSplitSize * 2); tree = fThreads[0]->GetAnalysisTree(); if (tree != nullptr) { @@ -409,8 +416,12 @@ void TRestProcessRunner::RunProcess() { } fOutputDataFile->cd(); - if (fEventTree != nullptr) fEventTree->Write(nullptr, kOverwrite); - if (fAnalysisTree != nullptr) fAnalysisTree->Write(nullptr, kOverwrite); + if (fEventTree != nullptr) { + fEventTree->Write(nullptr, kOverwrite); + } + if (fAnalysisTree != nullptr) { + fAnalysisTree->Write(nullptr, kOverwrite); + } // reset runner this->ResetRunTimes(); @@ -481,7 +492,7 @@ void TRestProcessRunner::RunProcess() { RESTEssential << "Waiting for processes to finish ..." << RESTendl; - while (1) { + while (true) { usleep(100000); bool finish = fThreads[0]->Finished(); for (int i = 1; i < fThreadNumber; i++) { @@ -557,8 +568,7 @@ void TRestProcessRunner::PauseMenu() { RESTLog << RESTendl; int menuupper = 15; int infobar = 11; - while (1) { - // Console::CursorUp(1); + while (true) { Console::ClearCurrentLine(); int b = Console::ReadKey(); // no need to press enter @@ -578,7 +588,7 @@ void TRestProcessRunner::PauseMenu() { RESTLog << "-" << RESTendl; RESTLog << RESTendl; RESTLog << RESTendl; - while (1) { + while (true) { Console::CursorUp(1); int c = Console::Read(); if (c != '\n') @@ -758,7 +768,7 @@ Int_t TRestProcessRunner::GetNextevtFunc(TRestEvent* targetevt, TRestAnalysisTre if (fProcessedEvents >= fEventsToProcess || targetevt == nullptr || fProcStatus == kStopping) { n = -1; } else { - if (fInputAnalysisStorage == false) { + if (!fInputAnalysisStorage) { n = fRunInfo->GetNextEvent(targetevt, nullptr); } else { n = fRunInfo->GetNextEvent(targetevt, targettree); @@ -782,8 +792,8 @@ Int_t TRestProcessRunner::GetNextevtFunc(TRestEvent* targetevt, TRestAnalysisTre void TRestProcessRunner::FillThreadEventFunc(TRestThread* t) { if (fSortOutputEvents) { // Make sure the thread has the minimum event id in the all the - // threads. Otherwise just wait. - while (1) { + // threads. Otherwise, just wait. + while (true) { bool smallest = true; for (TRestThread* otherT : fThreads) { if (otherT->Finished()) { @@ -848,8 +858,6 @@ void TRestProcessRunner::FillThreadEventFunc(TRestThread* t) { branchL->SetAddress(branchT->GetAddress()); } fEventTree->Fill(); - - // cout << t->GetOutputEvent()->GetID() << endl; } fProcessedEvents++; @@ -913,7 +921,7 @@ void TRestProcessRunner::FillThreadEventFunc(TRestThread* t) { } } - fOutputDataFile->Write(0, TObject::kOverwrite); + fOutputDataFile->Write(nullptr, TObject::kOverwrite); fOutputDataFile->Close(); delete fOutputDataFile; fOutputDataFile = newfile; @@ -947,8 +955,12 @@ void TRestProcessRunner::ConfigOutputFile() { #endif // write the last tree fOutputDataFile->cd(); - if (fEventTree != nullptr) fEventTree->Write(nullptr, kOverwrite); - if (fAnalysisTree != nullptr) fAnalysisTree->Write(nullptr, kOverwrite); + if (fEventTree != nullptr) { + fEventTree->Write(nullptr, kOverwrite); + } + if (fAnalysisTree != nullptr) { + fAnalysisTree->Write(nullptr, kOverwrite); + } // close file fOutputDataFile->Write(); @@ -964,30 +976,25 @@ void TRestProcessRunner::ConfigOutputFile() { // write metadata WriteProcessesMetadata(); + + // Write geometry + if (gGeoManager != nullptr) { + gGeoManager->Write("Geometry", TObject::kOverwrite); + } } /////////////////////////////////////////////// /// \brief Write process metadata to fOutputDataFile /// void TRestProcessRunner::WriteProcessesMetadata() { - // if (fRunInfo->GetInputFile() == nullptr) { - // if (!fRunInfo->GetOutputFile()) { - // // We are not ready yet to write - // return; - // } - // fRunInfo->cd(); - // } else - // fOutputDataFile->cd(); this->Write(nullptr, TObject::kWriteDelete); if (fRunInfo->GetFileProcess() != nullptr) { - // std::cout << "Run. Process-0. " << fRunInfo->GetFileProcess()->GetName() << std::endl; fRunInfo->GetFileProcess()->Write(nullptr, kOverwrite); } for (int i = 0; i < fProcessNumber; i++) { - // std::cout << "Thread. Process-" << i + 1 << fThreads[0]->GetProcess(i)->GetName() << std::endl; fThreads[0]->GetProcess(i)->Write(nullptr, kOverwrite); } } @@ -1006,7 +1013,7 @@ void TRestProcessRunner::MergeOutputFile() { for (int i = 0; i < fThreadNumber; i++) { TFile* f = fThreads[i]->GetOutputFile(); if (f != nullptr) { - f->Write(0, TObject::kOverwrite); + f->Write(nullptr, TObject::kOverwrite); f->Close(); } files_to_merge.push_back(f->GetName()); @@ -1017,8 +1024,6 @@ void TRestProcessRunner::MergeOutputFile() { } else { RESTError << "Output file: " << fOutputDataFileName << " is lost?" << RESTendl; } - // if (fRunInfo->GetInputFile() == nullptr) - // WriteProcessesMetadata(); } // tools @@ -1137,8 +1142,6 @@ void TRestProcessRunner::PrintProcessedEvents(Int_t rateE) { fflush(stdout); } - // CursorUp(4); - bytesAdded[poscalculated] = fRunInfo->GetBytesReaded() - bytesReaded_last; bytesReaded_last = fRunInfo->GetBytesReaded(); progAdded[poscalculated] = prog - prog_last; @@ -1169,16 +1172,8 @@ TRestEvent* TRestProcessRunner::GetInputEvent() { return fRunInfo->GetInputEvent TRestAnalysisTree* TRestProcessRunner::GetInputAnalysisTree() { return fRunInfo->GetAnalysisTree(); } void TRestProcessRunner::PrintMetadata() { - // cout.precision(10); TRestMetadata::PrintMetadata(); - /* - TRestStringOutput cout; - cout.setborder("||"); - cout.setorientation(1); - cout.setlength(100); - */ - string status; if (fProcStatus == kNormal) status = "Normal"; @@ -1188,16 +1183,12 @@ void TRestProcessRunner::PrintMetadata() { status = "Unknown"; RESTMetadata << "Status : " << status << RESTendl; - RESTMetadata << "Processesed events : " << fProcessedEvents << RESTendl; + RESTMetadata << "Processed events : " << fProcessedEvents << RESTendl; RESTMetadata << "Analysis tree branches : " << fNBranches << RESTendl; RESTMetadata << "Thread number : " << fThreadNumber << RESTendl; RESTMetadata << "Processes in each thread : " << fProcessNumber << RESTendl; RESTMetadata << "File auto split size: " << fFileSplitSize << RESTendl; RESTMetadata << "File compression level: " << fFileCompression << RESTendl; - // cout << "Input filename : " << fInputFilename << endl; - // cout << "Output filename : " << fOutputFilename << endl; - // cout << "Number of initial events : " << GetNumberOfEvents() << endl; - // cout << "Number of processed events : " << fProcessedEvents << endl; RESTMetadata << "******************************************" << RESTendl; RESTMetadata << RESTendl; RESTMetadata << RESTendl; diff --git a/source/framework/core/src/TRestRun.cxx b/source/framework/core/src/TRestRun.cxx index d08dba6b2..31f4da4e4 100644 --- a/source/framework/core/src/TRestRun.cxx +++ b/source/framework/core/src/TRestRun.cxx @@ -686,7 +686,7 @@ void TRestRun::ReadFileInfo(const string& filename) { int pos = -1; int pos1 = 0; int pos2 = 0; - while (1) { + while (true) { pos1 = format.find("[", pos + 1); pos2 = format.find("]", pos1); if (pos1 == -1 || pos2 == -1) { diff --git a/source/framework/core/src/TRestSystemOfUnits.cxx b/source/framework/core/src/TRestSystemOfUnits.cxx index a0fbc3e6a..ef01c59fc 100644 --- a/source/framework/core/src/TRestSystemOfUnits.cxx +++ b/source/framework/core/src/TRestSystemOfUnits.cxx @@ -310,8 +310,16 @@ TRestSystemOfUnits::TRestSystemOfUnits(string unitsStr) { } else { if (pos == unitsStr.size() - 1) { - RESTWarning << "last character inside \"" << unitsStr << "\" \"" << unitsStr[pos] + RESTWarning << "Last character inside \"" << unitsStr << "\" \"" << unitsStr[pos] << "\" unrecognized in unit definition!" << RESTendl; + + std::string lastChar = unitsStr.substr(pos, 1); + + if (isANumber(lastChar)) { + std::string tmpStr = unitsStr; + tmpStr.insert(pos, "^"); + RESTWarning << "Perhaps you meant: " << tmpStr << RESTendl; + } } pos++; diff --git a/source/framework/core/src/startup.cpp b/source/framework/core/src/startup.cpp index 78c19af1e..8e2481208 100644 --- a/source/framework/core/src/startup.cpp +++ b/source/framework/core/src/startup.cpp @@ -16,9 +16,9 @@ #include "TRestStringOutput.h" #include "TRestSystemOfUnits.h" #include "TRestTools.h" -#include "TRestVersion.h" + ////////////////////////////////////////////////////////////////////////// -/// This script initializes REST global variables in sequence to clearify +/// This script initializes REST global variables in sequence to clarify /// their dependency, therefore avoiding seg.fault during startup. All /// global variables in libRestTools, if depend on other global variable, /// should be placed here for initialization. @@ -77,11 +77,17 @@ struct __REST_CONST_INIT { char* _REST_PATH = getenv("REST_PATH"); char* _REST_USER = getenv("USER"); - char* _REST_USERHOME = getenv("HOME"); - // char* _REST_PATH = 0; - // char* _REST_USER = 0; - // char* _REST_USERHOME = 0; + // Should first look into REST_HOME, then HOME + char* _REST_USERHOME = getenv("REST_HOME"); + + if (_REST_USERHOME == nullptr) { + _REST_USERHOME = getenv("HOME"); + } else { + cout << "REST_HOME is set to " << _REST_USERHOME << endl; + // create the directory if it does not exist + std::filesystem::create_directories(_REST_USERHOME); + } #ifdef WIN32 if (_REST_PATH == nullptr) { @@ -272,6 +278,9 @@ string VectorToString(vector vec) { template vector StringToVector(string vec) { vector result; + + if (vec.empty()) return result; + if (vec[0] == '{' && vec[vec.size() - 1] == '}') { vec.erase(vec.begin()); vec.erase(vec.end() - 1); @@ -288,7 +297,9 @@ vector StringToVector(string vec) { } } else { - cout << "illegal format!" << endl; + cout << "Startup. StringToVector. Illegal format!" << endl; + cout << "The vector string is : " << vec << endl; + cout << "A vector should be defined using brackets and comma separated elements: {a,b,c,d}" << endl; return vector{}; } diff --git a/source/framework/tools/inc/TRestReflector.h b/source/framework/tools/inc/TRestReflector.h index 2a591175e..8607d52a9 100644 --- a/source/framework/tools/inc/TRestReflector.h +++ b/source/framework/tools/inc/TRestReflector.h @@ -306,9 +306,9 @@ class TRestReflector { /// Pointer to the corresponding TDataType helper, if the wrapped object is in data type bool is_data_type = false; /// If this object type wrapper is invalid - bool IsZombie(); + bool IsZombie() const; /// Deep copy the content of the wrapped object to `to`. - void operator>>(TRestReflector to); + void operator>>(const TRestReflector& to); /// Output overload by calling ToString(); friend std::ostream& operator<<(std::ostream& cin, TRestReflector ptr) { return cin << ptr.ToString(); } /// Get the value of the wrapped type, not recommended to use @@ -333,25 +333,25 @@ class TRestReflector { if (address != nullptr) *((T*)(address)) = val; } /// Convert the wrapped object to std::string - std::string ToString(); + std::string ToString() const; /// Set the value of the wrapped object from std::string - void ParseString(std::string str); + void ParseString(const std::string& str) const; /// Assembly a new object, and save its address. The old object will be destroied if not null void Assembly(); /// Destroy the current object. It will make the class to be zombie. - void Destroy(); + void Destroy() const; /// Print the Hex memory std::map of the wrappered object void PrintMemory(int bytepreline = 16); /// Find the class's datamember as TRestReflector object, including those from base class - TRestReflector GetDataMember(std::string name); + TRestReflector GetDataMember(const std::string& name); /// Get the i-th datamember of the class, ignoring those from base class TRestReflector GetDataMember(int ID); /// Get a list of the class's datamembers as a std::vector of std::string, including those from base class - std::vector GetListOfDataMembers(); + std::vector GetListOfDataMembers() const; /// Get the value of datamember as std::string. - std::string GetDataMemberValueString(std::string name); + std::string GetDataMemberValueString(const std::string& name); /// Get the number of data members of a class - int GetNumberOfDataMembers(); + int GetNumberOfDataMembers() const; /// Type conversion operator. With this, one can implicitly convert TRestReflector object to /// pointer of certain type. For example, `TRestEvent* eve = @@ -403,19 +403,19 @@ class TRestReflector { /////////////////////////////////////////////// /// \brief Assembly an object of type: typeName, returning the allocated memory address and size /// -TRestReflector Assembly(std::string typeName); +TRestReflector Assembly(const std::string& typeName); /////////////////////////////////////////////// /// \brief Wrap information an object of type: typeName, memory is not allocated /// -TRestReflector WrapType(std::string typeName); +TRestReflector WrapType(const std::string& typeName); /////////////////////////////////////////////// /// \brief Deep copy the content of object `from` to `to` /// /// Calls RESTVirtualConverter::CloneObj(). The actual methods are registered in converter.cpp /// If not registered, you can add it manually with AddConverter() macro -void CloneAny(TRestReflector from, TRestReflector to); +void CloneAny(const TRestReflector& from, const TRestReflector& to); }; // namespace REST_Reflection typedef REST_Reflection::TRestReflector RESTValue; diff --git a/source/framework/tools/inc/TRestStringHelper.h b/source/framework/tools/inc/TRestStringHelper.h index e9d804c02..879e5ee43 100644 --- a/source/framework/tools/inc/TRestStringHelper.h +++ b/source/framework/tools/inc/TRestStringHelper.h @@ -35,7 +35,7 @@ Float_t StringToFloat(std::string in); Double_t StringToDouble(std::string in); Int_t StringToInteger(std::string in); std::string IntegerToString(Int_t n, std::string format = "%d"); -std::string DoubleToString(Double_t d, std::string format = "%4.2lf"); +std::string DoubleToString(Double_t d, std::string format = "%8.6e"); Bool_t StringToBool(std::string booleanString); Long64_t StringToLong(std::string in); TVector3 StringTo3DVector(std::string in); diff --git a/source/framework/tools/inc/TRestTools.h b/source/framework/tools/inc/TRestTools.h index df1cf4360..58836b4bf 100644 --- a/source/framework/tools/inc/TRestTools.h +++ b/source/framework/tools/inc/TRestTools.h @@ -125,7 +125,7 @@ class TRestTools { static std::string Execute(std::string cmd); - static std::string DownloadRemoteFile(std::string remoteFile); + static std::string DownloadRemoteFile(const std::string& remoteFile); static int DownloadRemoteFile(std::string remoteFile, std::string localFile); static int UploadToServer(std::string localFile, std::string remoteFile, std::string methodUrl = ""); diff --git a/source/framework/tools/src/TRestReflector.cxx b/source/framework/tools/src/TRestReflector.cxx index fc34484d0..fccee751b 100644 --- a/source/framework/tools/src/TRestReflector.cxx +++ b/source/framework/tools/src/TRestReflector.cxx @@ -64,10 +64,10 @@ TRestReflector::TRestReflector(void* _address, const string& _type) { return; } - typeinfo = cl == 0 ? dt.typeinfo : cl->GetTypeInfo(); + typeinfo = cl == nullptr ? dt.typeinfo : cl->GetTypeInfo(); is_data_type = dt.size > 0; - size = cl == 0 ? dt.size : cl->Size(); - type = cl == 0 ? dt.name : cl->GetName(); + size = cl == nullptr ? dt.size : cl->Size(); + type = cl == nullptr ? dt.name : cl->GetName(); InitDictionary(); } @@ -87,9 +87,11 @@ void TRestReflector::Assembly() { } } -void TRestReflector::Destroy() { - if (address == nullptr) return; - if (onheap == false) { +void TRestReflector::Destroy() const { + if (address == nullptr) { + return; + } + if (!onheap) { // It can only delete/free objects on heap memory cout << "In TRestReflector::Destroy() : cannot free on stack memory!" << endl; return; @@ -119,11 +121,15 @@ void TRestReflector::PrintMemory(int bytepreline) { } } -void TRestReflector::operator>>(TRestReflector to) { CloneAny(*this, to); } +void TRestReflector::operator>>(const TRestReflector& to) { CloneAny(*this, to); } -string TRestReflector::ToString() { - if (type == "string") return *(string*)(address); - if (address == nullptr) return "null"; +string TRestReflector::ToString() const { + if (type == "string") { + return *(string*)(address); + } + if (address == nullptr) { + return "null"; + } RESTVirtualConverter* converter = RESTConverterMethodBase[typeinfo->hash_code()]; if (converter != nullptr) { return converter->ToString(address); @@ -132,7 +138,7 @@ string TRestReflector::ToString() { } } -void TRestReflector::ParseString(string str) { +void TRestReflector::ParseString(const string& str) const { if (type == "string") { *(string*)(address) = str; } else { @@ -146,7 +152,9 @@ void TRestReflector::ParseString(string str) { } int TRestReflector::InitDictionary() { - if (is_data_type) return 0; + if (is_data_type) { + return 0; + } if (cl != nullptr) { if (cl->GetCollectionProxy() && dynamic_cast(cl->GetCollectionProxy())) { @@ -159,12 +167,12 @@ int TRestReflector::InitDictionary() { } } - if (type == "" || size == 0 || cl == nullptr) { + if (type.empty() || size == 0 || cl == nullptr) { cout << "Error in CreateDictionary: object is zombie!" << endl; return -1; } - int pos = ((string)type).find("<"); + int pos = ((string)type).find('<'); string basetype = ((string)type).substr(0, pos); vector stltypes{"vector", "list", "map", "set", "array", "deque"}; @@ -190,7 +198,6 @@ int TRestReflector::InitDictionary() { if (TRestTools::fileExists(sofilename)) { cout << "Loading external dictionary for: \"" << type << "\":" << endl; cout << sofilename << endl; - } else { // we create a new library of dictionary for that type if (!TRestTools::isPathWritable(REST_USER_PATH)) { @@ -201,8 +208,7 @@ int TRestReflector::InitDictionary() { << endl; return -1; } - int z = system(Form("mkdir -p %s/AddonDict", REST_USER_PATH.c_str())); - if (z != 0) { + if (system(Form("mkdir -p %s/AddonDict", REST_USER_PATH.c_str())) != 0) { cout << "mkdir failed to create directory" << endl; return -1; } @@ -225,26 +231,15 @@ int TRestReflector::InitDictionary() { cout << "Creating external dictionary for: \"" << type << "\":" << endl; cout << sofilename << endl; - // cout << Form("rootcling -f %s -c %s", outfilename.c_str(), infilename.c_str()) << endl; - int a = system(Form("rootcling -f %s -c %s", cxxfilename.c_str(), linkdeffilename.c_str())); - if (a != 0) { + if (system(Form("rootcling -f %s -c %s", cxxfilename.c_str(), linkdeffilename.c_str())) != 0) { cout << "rootcling failed to generate dictionary" << endl; return -1; } - int b = - system(Form("gcc %s `root-config --cflags` " + if (system(Form("gcc %s `root-config --cflags` " "`root-config --libs` -lGui -lGeom -lGdml -lMinuit -L/usr/lib64 " "-lstdc++ -shared -fPIC -o %s", - cxxfilename.c_str(), sofilename.c_str())); - - // int c = - // system(Form("gcc %s/lib/AddonDict/*.cxx -std=c++11 -I`root-config --incdir` " - // "`root-config --libs` -lGui -lEve -lGeom -lMathMore -lGdml -lMinuit -L/usr/lib64 " - // "-lstdc++ -shared -fPIC -o %s/lib/AddonDict/libRestAddonDict.so", - // restpath, restpath)); - - if (b != 0 /*|| c != 0*/) { + cxxfilename.c_str(), sofilename.c_str())) != 0) { cout << "gcc failed to generate library for the dictionary" << endl; return -1; } @@ -258,19 +253,19 @@ int TRestReflector::InitDictionary() { return 0; } -bool TRestReflector::IsZombie() { - return (type == "" || address == 0 || size == 0 || (cl == 0 && !is_data_type)); +bool TRestReflector::IsZombie() const { + return (type.empty() || address == nullptr || size == 0 || (cl == nullptr && !is_data_type)); } -TRestReflector Assembly(string typeName) { +TRestReflector Assembly(const string& typeName) { TRestReflector ptr = WrapType(typeName); ptr.Assembly(); return ptr; } -TRestReflector WrapType(string type) { return TRestReflector(0, type); } +TRestReflector WrapType(const string& type) { return TRestReflector(nullptr, type); } -void CloneAny(TRestReflector from, TRestReflector to) { +void CloneAny(const TRestReflector& from, const TRestReflector& to) { if (from.IsZombie() || to.IsZombie()) { cout << "In TRestReflector::CloneTo() : the ptr is zombie! " << endl; return; @@ -288,24 +283,9 @@ void CloneAny(TRestReflector from, TRestReflector to) { } else { cout << "Method for cloning type: \"" << from.type << "\" has not been registered!" << endl; } - // if (from.cl == nullptr) { - // memcpy(to.address, from.address, from.size); - //} else { - // TBufferFile buffer(TBuffer::kWrite); - - // buffer.MapObject(from.address, from.cl); // register obj in map to handle self reference - // from.cl->Streamer(from.address, buffer); - - // buffer.SetReadMode(); - // buffer.ResetMap(); - // buffer.SetBufferOffset(0); - - // buffer.MapObject(to.address, to.cl); // register obj in map to handle self reference - // to.cl->Streamer(to.address, buffer); - //} } -TRestReflector TRestReflector::GetDataMember(string name) { +TRestReflector TRestReflector::GetDataMember(const string& name) { if (cl != nullptr) { TDataMember* mem = cl->GetDataMember(name.c_str()); if (mem == nullptr) { @@ -349,7 +329,7 @@ TRestReflector TRestReflector::GetDataMember(int ID) { return TRestReflector(); } -vector TRestReflector::GetListOfDataMembers() { +vector TRestReflector::GetListOfDataMembers() const { vector dataMembers; // add datamembers from base class first @@ -380,7 +360,7 @@ vector TRestReflector::GetListOfDataMembers() { return dataMembers; } -string TRestReflector::GetDataMemberValueString(string name) { +string TRestReflector::GetDataMemberValueString(const string& name) { TRestReflector member = GetDataMember(name); if (!member.IsZombie()) { return member.ToString(); @@ -388,7 +368,7 @@ string TRestReflector::GetDataMemberValueString(string name) { return ""; } -int TRestReflector::GetNumberOfDataMembers() { +int TRestReflector::GetNumberOfDataMembers() const { if (cl != nullptr) { return cl->GetNdata(); } diff --git a/source/framework/tools/src/TRestTools.cxx b/source/framework/tools/src/TRestTools.cxx index aa94f34c8..0af3237fa 100644 --- a/source/framework/tools/src/TRestTools.cxx +++ b/source/framework/tools/src/TRestTools.cxx @@ -49,6 +49,8 @@ #include #include +#include + #ifdef USE_Curl #include #endif @@ -711,11 +713,9 @@ bool TRestTools::isDataSet(const std::string& filename) { /////////////////////////////////////////////// /// \brief Returns true if **filename** is an *http* address. /// -bool TRestTools::isURL(const string& filename) { - if (filename.find("http") == 0) { - return true; - } - return false; +bool TRestTools::isURL(const string& s) { + std::regex pattern("^https?://(.+)"); + return std::regex_match(s, pattern); } /////////////////////////////////////////////// @@ -1049,19 +1049,19 @@ std::istream& TRestTools::GetLine(std::istream& is, std::string& t) { /// will be used, including scp, wget. Downloads to REST_USER_PATH + "/download/" + filename /// by default. /// -std::string TRestTools::DownloadRemoteFile(string url) { - string purename = TRestTools::GetPureFileName(url); - if (purename == "") { - cout << "error! (TRestTools::DownloadRemoteFile): url is not a file!" << endl; - cout << "please specify a concrete file name in this url" << endl; - cout << "url: " << url << endl; +std::string TRestTools::DownloadRemoteFile(const string& url) { + string pureName = TRestTools::GetPureFileName(url); + if (pureName.empty()) { + RESTWarning << "error! (TRestTools::DownloadRemoteFile): url is not a file!" << RESTendl; + RESTWarning << "please specify a concrete file name in this url" << RESTendl; + RESTWarning << "url: " << url << RESTendl; return ""; } if (url.find("local:") == 0) { return Replace(url, "local:", ""); } else { - string fullpath = REST_USER_PATH + "/download/" + purename; + string fullpath = REST_USER_PATH + "/download/" + pureName; int out; int attempts = 10; do { @@ -1070,14 +1070,12 @@ std::string TRestTools::DownloadRemoteFile(string url) { RESTWarning << "Retrying download in 5 seconds" << RESTendl; std::this_thread::sleep_for(std::chrono::seconds(5)); } else if (attempts < 10) { - RESTSuccess << "Download suceeded after " << 10 - attempts << " attempts" << RESTendl; + RESTSuccess << "Download succeeded after " << 10 - attempts << " attempts" << RESTendl; } attempts--; } while (out == 1024 && attempts > 0); - if (out == 0) { - return fullpath; - } else if (TRestTools::fileExists(fullpath)) { + if (out == 0 || TRestTools::fileExists(fullpath)) { return fullpath; } else { return "";