From 05748ac7f002f2c15591d0a7d6f636685858d9ce Mon Sep 17 00:00:00 2001 From: lobis Date: Sun, 14 Feb 2021 11:55:26 +0100 Subject: [PATCH 1/4] added capture volumes --- ...eant4CosmicNeutronTaggingAnalysisProcess.h | 135 ++++++++ ...nt4CosmicNeutronTaggingAnalysisProcess.cxx | 296 ++++++++++++++++++ 2 files changed, 431 insertions(+) create mode 100644 inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h create mode 100644 src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx diff --git a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h new file mode 100644 index 00000000..8ebfbccb --- /dev/null +++ b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h @@ -0,0 +1,135 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +#ifndef RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess +#define RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess + +#include +#include + +#include "TRestEventProcess.h" + +class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess { + private: + /// A pointer to the specific TRestGeant4Event input + TRestGeant4Event* fInputG4Event; //! + /// A pointer to the specific TRestGeant4Event output + TRestGeant4Event* fOutputG4Event; //! + /// A pointer to the simulation metadata information accessible to TRestRun + TRestGeant4Metadata* fG4Metadata; //! + + std::vector fVetoVolumeIds; //! + std::vector fCaptureVolumeIds; //! + string fVetoKeyword = ""; //! + string fCaptureKeyword = ""; //! + std::vector fVetoGroupKeywords; //! + std::map> fVetoGroupVolumeNames; //! + std::vector fQuenchingFactors; //! + + void InitFromConfigFile(); + void Initialize(); + void LoadDefaultConfig(); + + // clean string (https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring) + inline std::string& rtrim(std::string& s, const char* t = " \t\n\r\f\v") { + s.erase(s.find_last_not_of(t) + 1); + return s; + } + // trim from beginning of string (left) + inline std::string& ltrim(std::string& s, const char* t = " \t\n\r\f\v") { + s.erase(0, s.find_first_not_of(t)); + return s; + } + // trim from both ends of string (right then left) + inline std::string& trim(std::string& s, const char* t = " \t\n\r\f\v") { return ltrim(rtrim(s, t), t); } + + // final clean string: trim and UPPER + inline std::string& clean_string(std::string& s) { + s = trim(s); + std::transform(s.begin(), s.end(), s.begin(), ::tolower); + return s; + } + + protected: + // add here the members of your event process + + public: + any GetInputEvent() { return fInputG4Event; } + any GetOutputEvent() { return fOutputG4Event; } + + void InitProcess(); + TRestEvent* ProcessEvent(TRestEvent* eventInput); + void EndProcess(); + + void LoadConfig(std::string cfgFilename, std::string name = ""); + + /// It prints out the process parameters stored in the metadata structure + void PrintMetadata() { + BeginPrintProcess(); + + debug << "VETO KEYWORD: " << fVetoKeyword << endl; + debug << endl; + + debug << "VETO GROUP KEYWORDS:" << endl; + for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { + debug << "\t" << fVetoGroupKeywords[i] << endl; + } + debug << endl; + + debug << "Found " << fVetoVolumeIds.size() << " veto volumes:" << endl; + for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fVetoVolumeIds[i]) << endl; + } + debug << endl; + + // capture volumes + + debug << "CAPTURE KEYWORD: " << fCaptureKeyword << endl; + debug << endl; + + debug << "Found " << fCaptureVolumeIds.size() << " Capture volumes:" << endl; + for (unsigned int i = 0; i < fCaptureVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fCaptureVolumeIds[i]) << endl; + } + debug << endl; + + debug << "QUENCHING FACTORS (" << fQuenchingFactors.size() << " Total)" << endl; + for (unsigned int i = 0; i < fQuenchingFactors.size(); i++) { + debug << "\t" << fQuenchingFactors[i] << endl; + } + debug << endl; + + EndPrintProcess(); + } + + /// Returns a new instance of this class + TRestEventProcess* Maker() { return new TRestGeant4CosmicNeutronTaggingAnalysisProcess; } + /// Returns the name of this process + TString GetProcessName() { return (TString) "geant4VetoAnalysis"; } + + TRestGeant4CosmicNeutronTaggingAnalysisProcess(); + TRestGeant4CosmicNeutronTaggingAnalysisProcess(char* cfgFileName); + ~TRestGeant4CosmicNeutronTaggingAnalysisProcess(); + + ClassDef(TRestGeant4CosmicNeutronTaggingAnalysisProcess, 1); +}; +#endif // RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess diff --git a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx new file mode 100644 index 00000000..bf83e62e --- /dev/null +++ b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx @@ -0,0 +1,296 @@ +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// TODO: fill documentation +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2021-February: Implementation. +/// +/// \class TRestGeant4CosmicNeutronTaggingAnalysisProcess +/// \author Luis Obis +/// +///
+/// + +#include "TRestGeant4CosmicNeutronTaggingAnalysisProcess.h" +using namespace std; + +ClassImp(TRestGeant4CosmicNeutronTaggingAnalysisProcess); + +TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess() { + Initialize(); +} + +TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess( + char* cfgFileName) { + Initialize(); + if (LoadConfigFromFile(cfgFileName)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestGeant4CosmicNeutronTaggingAnalysisProcess::~TRestGeant4CosmicNeutronTaggingAnalysisProcess() { + delete fOutputG4Event; +} + +/////////////////////////////////////////////// +/// \brief Function to load the default config in absence of RML input +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); } + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define the +/// section name +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Initialize() { + fG4Metadata = NULL; + + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fInputG4Event = NULL; + fOutputG4Event = new TRestGeant4Event(); +} + +/////////////////////////////////////////////// +/// \brief Function to load the configuration from an external configuration +/// file. +/// +/// If no configuration path is defined in TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// correspondig TRestGeant4CosmicNeutronTaggingAnalysisProcess section inside the RML. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadConfig(std::string cfgFilename, std::string name) { + if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Process initialization. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { + fG4Metadata = GetMetadata(); + + // CAREFUL THIS METHOD IS CALLED TWICE! + // we need to reset these variables to zero + fVetoVolumeIds.clear(); + fVetoGroupVolumeNames.clear(); + fCaptureVolumeIds.clear(); + + // get "veto" volumes + for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); + volume_name = clean_string(volume_name); + if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { + fVetoVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { + fCaptureVolumeIds.push_back(i); + } + } + + // veto groups (fill fVetoGroupVolumeNames) + for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { + string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); + fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; + for (int& id : fVetoVolumeIds) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + volume_name = clean_string(volume_name); + if (volume_name.find(veto_group_keyword) != string::npos) { + fVetoGroupVolumeNames[veto_group_keyword].push_back( + (string)fG4Metadata->GetActiveVolumeName(id)); + } + } + } + + PrintMetadata(); +} + +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEvent* evInput) { + fInputG4Event = (TRestGeant4Event*)evInput; + *fOutputG4Event = *((TRestGeant4Event*)evInput); + + std::map volume_energy_map; + + for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { + int id = fVetoVolumeIds[i]; + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + + Double_t energy = fOutputG4Event->GetEnergyDepositedInVolume(id); + volume_energy_map[volume_name] = energy; + } + + Double_t energy_veto_max = 0; + for (const auto& pair : volume_energy_map) { + Double_t veto_energy = pair.second; + SetObservableValue(pair.first + "VolumeEDep", veto_energy); + if (veto_energy > energy_veto_max) { + energy_veto_max = veto_energy; + }; + } + SetObservableValue("vetoAllEVetoMax", energy_veto_max); + + // veto groups + for (const auto& pair : fVetoGroupVolumeNames) { + Double_t energy_veto_max_group = 0; + for (unsigned int i = 0; i < pair.second.size(); i++) { + string volume_name = pair.second[i]; + Double_t veto_energy = volume_energy_map[volume_name]; + if (veto_energy > energy_veto_max_group) { + energy_veto_max_group = veto_energy; + }; + } + // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) + string group_name; + for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { + if (it == pair.first.cbegin()) { + group_name += std::toupper(*it); + } else { + group_name += std::tolower(*it); + } + } + SetObservableValue("vetoGroup" + group_name + "EVetoMax", energy_veto_max_group); + } + + // + + for (const auto& quenching_factor : fQuenchingFactors) { + string quenching_factor_string = std::to_string(quenching_factor); + // replace "." in string by "_" because its gives very strange problems + quenching_factor_string = + quenching_factor_string.replace(quenching_factor_string.find("."), sizeof(".") - 1, "_"); + volume_energy_map.clear(); + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + for (const auto& id : fVetoVolumeIds) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + + if (particle_name == "e-" || particle_name == "e+" || particle_name == "gamma") { + // no quenching factor + volume_energy_map[volume_name] += track->GetEnergyInVolume(id); + } else { + // apply quenching factor + volume_energy_map[volume_name] += quenching_factor * track->GetEnergyInVolume(id); + } + } + } + + Double_t energy_veto_max = 0; + for (const auto& pair : volume_energy_map) { + Double_t veto_energy = pair.second; + SetObservableValue(pair.first + "VolumeEDep" + "Qf" + quenching_factor_string, veto_energy); + if (veto_energy > energy_veto_max) { + energy_veto_max = veto_energy; + }; + } + SetObservableValue(string("vetoAllEVetoMax") + "Qf" + quenching_factor_string, energy_veto_max); + + // veto groups + for (const auto& pair : fVetoGroupVolumeNames) { + Double_t energy_veto_max_group = 0; + for (unsigned int i = 0; i < pair.second.size(); i++) { + string volume_name = pair.second[i]; + Double_t veto_energy = volume_energy_map[volume_name]; + if (veto_energy > energy_veto_max_group) { + energy_veto_max_group = veto_energy; + }; + } + // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) + string group_name; + for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { + if (it == pair.first.cbegin()) { + group_name += std::toupper(*it); + } else { + group_name += std::tolower(*it); + } + } + SetObservableValue("vetoGroup" + group_name + "EVetoMax" + "Qf" + quenching_factor_string, + energy_veto_max_group); + } + } + + return fOutputG4Event; +} + +/////////////////////////////////////////////// +/// \brief Function to include required actions after all events have been +/// processed. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::EndProcess() { + // Function to be executed once at the end of the process + // (after all events have been processed) + + // Start by calling the EndProcess function of the abstract class. + // Comment this if you don't want it. + // TRestEventProcess::EndProcess(); +} + +/////////////////////////////////////////////// +/// \brief Function to read input parameters from the RML +/// TRestGeant4CosmicNeutronTaggingAnalysisProcess metadata section +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitFromConfigFile() { + // word to identify active volume as veto (default = "veto" e.g. "vetoTop") + string veto_keyword = GetParameter("vetoKeyword", "veto"); + fVetoKeyword = clean_string(veto_keyword); + // comma separated tags: "top, bottom, ..." + string veto_group_keywords = GetParameter("vetoGroupKeywords", ""); + stringstream ss(veto_group_keywords); + while (ss.good()) { + string substr; + getline(ss, substr, ','); + fVetoGroupKeywords.push_back(clean_string(substr)); + } + + // word to identify active volume as capture sheet (cadmium, default = "sheet" e.g. + // "scintillatorSheetTop1of4") + string capture_keyword = GetParameter("captureKeyword", "sheet"); + fCaptureKeyword = clean_string(capture_keyword); + + // comma separated quenching factors: "0.15, 1.00, ..." + string quenching_factors = GetParameter("vetoQuenchingFactors", "-1"); + stringstream ss_qf(quenching_factors); + while (ss_qf.good()) { + string substr; + getline(ss_qf, substr, ','); + substr = clean_string(substr); + Float_t quenching_factor = (Float_t)std::atof(substr.c_str()); + if (quenching_factor > 1 || quenching_factor < 0) { + cout << "ERROR: quenching factor must be between 0 and 1" << endl; + continue; + } + fQuenchingFactors.push_back(quenching_factor); + } +} From bd08013883ed3727270818c81e2b8b7f2d27a6e9 Mon Sep 17 00:00:00 2001 From: lobis Date: Sun, 14 Feb 2021 15:05:53 +0100 Subject: [PATCH 2/4] added neutron analysis process --- ...eant4CosmicNeutronTaggingAnalysisProcess.h | 46 ++++- ...nt4CosmicNeutronTaggingAnalysisProcess.cxx | 186 +++++++++++++++++- 2 files changed, 225 insertions(+), 7 deletions(-) diff --git a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h index 8ebfbccb..0a106a7d 100644 --- a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h +++ b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h @@ -37,18 +37,47 @@ class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess /// A pointer to the simulation metadata information accessible to TRestRun TRestGeant4Metadata* fG4Metadata; //! - std::vector fVetoVolumeIds; //! - std::vector fCaptureVolumeIds; //! + std::vector fVetoVolumeIds; //! + std::vector fCaptureVolumeIds; //! + std::vector fShieldingVolumeIds; //! + string fVetoKeyword = ""; //! string fCaptureKeyword = ""; //! + string fShieldingKeyword = ""; //! std::vector fVetoGroupKeywords; //! std::map> fVetoGroupVolumeNames; //! std::vector fQuenchingFactors; //! + // neutrons that undergo neutron capture + Int_t fNeutronsCapturedNumber; //! + std::vector fNeutronsCapturedPosX; //! + std::vector fNeutronsCapturedPosY; //! + std::vector fNeutronsCapturedPosZ; //! + std::vector fNeutronsCapturedIsCaptureVolume; //! + std::vector fNeutronsCapturedProductionE; //! + + // gammas that are produced from neutron capture + Int_t fGammasNeutronCaptureNumber; //! + std::vector fGammasNeutronCapturePosX; //! + std::vector fGammasNeutronCapturePosY; //! + std::vector fGammasNeutronCapturePosZ; //! + std::vector fGammasNeutronCaptureIsCaptureVolume; //! + std::vector fGammasNeutronCaptureProductionE; //! + + // secondary neutrons that exit the lead shielding + Int_t fSecondaryNeutronsShieldingNumber; //! + std::vector fSecondaryNeutronsShieldingExitPosX; //! + std::vector fSecondaryNeutronsShieldingExitPosY; //! + std::vector fSecondaryNeutronsShieldingExitPosZ; //! + std::vector fSecondaryNeutronsShieldingIsCaptured; //! + std::vector fSecondaryNeutronsShieldingIsCapturedInCaptureVolume; //! + std::vector fSecondaryNeutronsShieldingProductionE; //! + std::vector fSecondaryNeutronsShieldingExitE; //! + void InitFromConfigFile(); void Initialize(); void LoadDefaultConfig(); - + void Reset(); // clean string (https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring) inline std::string& rtrim(std::string& s, const char* t = " \t\n\r\f\v") { s.erase(s.find_last_not_of(t) + 1); @@ -112,6 +141,17 @@ class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess } debug << endl; + // shielding volume/s + + debug << "SHIELDING KEYWORD: " << fShieldingKeyword << endl; + debug << endl; + + debug << "Found " << fShieldingVolumeIds.size() << " Shielding volumes:" << endl; + for (unsigned int i = 0; i < fShieldingVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fShieldingVolumeIds[i]) << endl; + } + debug << endl; + debug << "QUENCHING FACTORS (" << fQuenchingFactors.size() << " Total)" << endl; for (unsigned int i = 0; i < fQuenchingFactors.size(); i++) { debug << "\t" << fQuenchingFactors[i] << endl; diff --git a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx index bf83e62e..5c4667ad 100644 --- a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx +++ b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx @@ -102,10 +102,7 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { // CAREFUL THIS METHOD IS CALLED TWICE! // we need to reset these variables to zero - fVetoVolumeIds.clear(); - fVetoGroupVolumeNames.clear(); - fCaptureVolumeIds.clear(); - + Reset(); // get "veto" volumes for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); @@ -114,6 +111,8 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { fVetoVolumeIds.push_back(i); } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { fCaptureVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fShieldingKeyword)) != string::npos) { + fShieldingVolumeIds.push_back(i); } } @@ -134,6 +133,34 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { PrintMetadata(); } +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Reset() { + fVetoVolumeIds.clear(); + fVetoGroupVolumeNames.clear(); + fCaptureVolumeIds.clear(); + + fNeutronsCapturedNumber = 0; + fNeutronsCapturedPosX.clear(); + fNeutronsCapturedPosY.clear(); + fNeutronsCapturedPosZ.clear(); + fNeutronsCapturedIsCaptureVolume.clear(); + fNeutronsCapturedProductionE.clear(); + + fGammasNeutronCaptureNumber = 0; + fGammasNeutronCapturePosX.clear(); + fGammasNeutronCapturePosY.clear(); + fGammasNeutronCapturePosZ.clear(); + fGammasNeutronCaptureIsCaptureVolume.clear(); + fGammasNeutronCaptureProductionE.clear(); + + fSecondaryNeutronsShieldingNumber = 0; + fSecondaryNeutronsShieldingExitPosX.clear(); + fSecondaryNeutronsShieldingExitPosY.clear(); + fSecondaryNeutronsShieldingExitPosZ.clear(); + fSecondaryNeutronsShieldingIsCaptured.clear(); + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.clear(); + fSecondaryNeutronsShieldingProductionE.clear(); + fSecondaryNeutronsShieldingExitE.clear(); +} /////////////////////////////////////////////// /// \brief The main processing event function /// @@ -141,6 +168,7 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv fInputG4Event = (TRestGeant4Event*)evInput; *fOutputG4Event = *((TRestGeant4Event*)evInput); + Reset(); std::map volume_energy_map; for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { @@ -241,6 +269,151 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv } } + std::set neutronsCaptured = {}; + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "neutron") { + auto hits = track->GetHits(); + for (int j = 0; j < hits->GetNumberOfHits(); j++) { + string process_name = (string)track->GetProcessName(hits->GetProcess(j)); + if (process_name == "nCapture") { + cout << "Neutron capture!!!!!! " << particle_name << "trackId " << track->GetTrackID() + << " hit " << j << endl; + // track->PrintTrack(); + hits->PrintHits(j + 1); + + neutronsCaptured.insert(track->GetTrackID()); + + fNeutronsCapturedNumber += 1; + fNeutronsCapturedPosX.push_back(hits->GetX(j)); + fNeutronsCapturedPosY.push_back(hits->GetY(j)); + fNeutronsCapturedPosZ.push_back(hits->GetZ(j)); + + Int_t volumeId = hits->GetVolumeId(j); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + fNeutronsCapturedIsCaptureVolume.push_back(isCaptureVolume); + fNeutronsCapturedProductionE.push_back(track->GetKineticEnergy()); + } + } + } + } + + SetObservableValue("neutronsCapturedNumber", fNeutronsCapturedNumber); + SetObservableValue("neutronsCapturedPosX", fNeutronsCapturedPosX); + SetObservableValue("neutronsCapturedPosY", fNeutronsCapturedPosY); + SetObservableValue("neutronsCapturedPosZ", fNeutronsCapturedPosZ); + SetObservableValue("neutronsCapturedIsCaptureVolume", fNeutronsCapturedIsCaptureVolume); + SetObservableValue("neutronsCapturedProductionE", fNeutronsCapturedProductionE); + + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "gamma") { + // check if gamma is child of captured neutron + Int_t parent = track->GetParentID(); + if (neutronsCaptured.count(parent) > 0) { + auto hits = track->GetHits(); + + fGammasNeutronCaptureNumber += 1; + fGammasNeutronCapturePosX.push_back(hits->GetX(0)); + fGammasNeutronCapturePosY.push_back(hits->GetY(0)); + fGammasNeutronCapturePosZ.push_back(hits->GetZ(0)); + + Int_t volumeId = hits->GetVolumeId(0); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + fGammasNeutronCaptureIsCaptureVolume.push_back(isCaptureVolume); + fGammasNeutronCaptureProductionE.push_back(track->GetKineticEnergy()); + + cout << "gamma capture" << endl; + + hits->PrintHits(1); + } + } + } + + SetObservableValue("gammasNeutronCaptureNumber", fGammasNeutronCaptureNumber); + SetObservableValue("gammasNeutronCapturePosX", fGammasNeutronCapturePosX); + SetObservableValue("gammasNeutronCapturePosY", fGammasNeutronCapturePosY); + SetObservableValue("gammasNeutronCapturePosZ", fGammasNeutronCapturePosZ); + SetObservableValue("gammasNeutronCaptureIsCaptureVolume", fGammasNeutronCaptureIsCaptureVolume); + SetObservableValue("gammasNeutronCaptureProductionE", fGammasNeutronCaptureProductionE); + + std::set secondaryNeutrons = {}; // avoid counting twice + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "neutron" && track->GetParentID() != 0) { // not consider primary + // check if neutron exits shielding + auto hits = track->GetHits(); + for (int j = 0; j < hits->GetNumberOfHits(); j++) { + string process_name = (string)track->GetProcessName(hits->GetProcess(j)); + if (process_name == "Transportation") { + for (const auto& id : fShieldingVolumeIds) { + if (hits->GetVolumeId(j) == id) { + // transportation and shielding == exits shielding + if (secondaryNeutrons.count(track->GetTrackID()) == 0) { + // first time adding this secondary neutron + secondaryNeutrons.insert(track->GetTrackID()); + } else { + continue; + } + fSecondaryNeutronsShieldingNumber += 1; + fSecondaryNeutronsShieldingExitPosX.push_back(hits->GetX(j)); + fSecondaryNeutronsShieldingExitPosY.push_back(hits->GetY(j)); + fSecondaryNeutronsShieldingExitPosZ.push_back(hits->GetZ(j)); + + Int_t volumeId = hits->GetVolumeId(j); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + Int_t isCaptured = 0; + if (neutronsCaptured.count(track->GetTrackID()) > 0) { + isCaptured = 1; + } + fSecondaryNeutronsShieldingIsCaptured.push_back(isCaptured); + if (isCaptured) + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back( + isCaptureVolume); + else { + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back(0); + } + + fSecondaryNeutronsShieldingProductionE.push_back(track->GetKineticEnergy()); + fSecondaryNeutronsShieldingExitE.push_back(hits->GetKineticEnergy(j)); + } + } + } + } + } + } + + SetObservableValue("fSecondaryNeutronsShieldingNumber", fSecondaryNeutronsShieldingNumber); + SetObservableValue("fSecondaryNeutronsShieldingExitPosX", fSecondaryNeutronsShieldingExitPosX); + SetObservableValue("fSecondaryNeutronsShieldingExitPosY", fSecondaryNeutronsShieldingExitPosY); + SetObservableValue("fSecondaryNeutronsShieldingExitPosZ", fSecondaryNeutronsShieldingExitPosZ); + SetObservableValue("fSecondaryNeutronsShieldingIsCaptured", fSecondaryNeutronsShieldingIsCaptured); + SetObservableValue("fSecondaryNeutronsShieldingIsCapturedInCaptureVolume", + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume); + SetObservableValue("fSecondaryNeutronsShieldingProductionE", fSecondaryNeutronsShieldingProductionE); + SetObservableValue("fSecondaryNeutronsShieldingExitE", fSecondaryNeutronsShieldingExitE); + return fOutputG4Event; } @@ -279,6 +452,11 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitFromConfigFile() { string capture_keyword = GetParameter("captureKeyword", "sheet"); fCaptureKeyword = clean_string(capture_keyword); + // word to identify active volume as shielding + + string shielding_keyword = GetParameter("shieldingKeyword", "shielding"); + fShieldingKeyword = clean_string(shielding_keyword); + // comma separated quenching factors: "0.15, 1.00, ..." string quenching_factors = GetParameter("vetoQuenchingFactors", "-1"); stringstream ss_qf(quenching_factors); From 37341f340d67ec5c1a3f8b57083d933e08e8b216 Mon Sep 17 00:00:00 2001 From: lobis Date: Sun, 14 Feb 2021 21:17:39 +0100 Subject: [PATCH 3/4] added more neutron observables --- ...eant4CosmicNeutronTaggingAnalysisProcess.h | 18 ++- ...nt4CosmicNeutronTaggingAnalysisProcess.cxx | 136 +++++++++++++----- 2 files changed, 113 insertions(+), 41 deletions(-) diff --git a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h index 0a106a7d..9791f0f2 100644 --- a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h +++ b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h @@ -49,12 +49,18 @@ class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess std::vector fQuenchingFactors; //! // neutrons that undergo neutron capture - Int_t fNeutronsCapturedNumber; //! - std::vector fNeutronsCapturedPosX; //! - std::vector fNeutronsCapturedPosY; //! - std::vector fNeutronsCapturedPosZ; //! - std::vector fNeutronsCapturedIsCaptureVolume; //! - std::vector fNeutronsCapturedProductionE; //! + Int_t fNeutronsCapturedNumber; //! + std::vector fNeutronsCapturedPosX; //! + std::vector fNeutronsCapturedPosY; //! + std::vector fNeutronsCapturedPosZ; //! + std::vector fNeutronsCapturedIsCaptureVolume; //! + std::vector fNeutronsCapturedProductionE; //! + std::vector fNeutronsCapturedEDepByNeutron; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildren; //! + std::vector fNeutronsCapturedEDepByNeutronInVeto; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVeto; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin; //! // gammas that are produced from neutron capture Int_t fGammasNeutronCaptureNumber; //! diff --git a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx index 5c4667ad..01c7a9ca 100644 --- a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx +++ b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx @@ -104,28 +104,30 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { // we need to reset these variables to zero Reset(); // get "veto" volumes - for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); - volume_name = clean_string(volume_name); - if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { - fVetoVolumeIds.push_back(i); - } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { - fCaptureVolumeIds.push_back(i); - } else if (volume_name.find(clean_string(fShieldingKeyword)) != string::npos) { - fShieldingVolumeIds.push_back(i); + if (fVetoVolumeIds.empty()) { + for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); + volume_name = clean_string(volume_name); + if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { + fVetoVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { + fCaptureVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fShieldingKeyword)) != string::npos) { + fShieldingVolumeIds.push_back(i); + } } - } - // veto groups (fill fVetoGroupVolumeNames) - for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { - string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); - fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; - for (int& id : fVetoVolumeIds) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); - volume_name = clean_string(volume_name); - if (volume_name.find(veto_group_keyword) != string::npos) { - fVetoGroupVolumeNames[veto_group_keyword].push_back( - (string)fG4Metadata->GetActiveVolumeName(id)); + // veto groups (fill fVetoGroupVolumeNames) + for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { + string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); + fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; + for (int& id : fVetoVolumeIds) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + volume_name = clean_string(volume_name); + if (volume_name.find(veto_group_keyword) != string::npos) { + fVetoGroupVolumeNames[veto_group_keyword].push_back( + (string)fG4Metadata->GetActiveVolumeName(id)); + } } } } @@ -134,16 +136,23 @@ void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { } void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Reset() { + /* fVetoVolumeIds.clear(); fVetoGroupVolumeNames.clear(); fCaptureVolumeIds.clear(); - + */ fNeutronsCapturedNumber = 0; fNeutronsCapturedPosX.clear(); fNeutronsCapturedPosY.clear(); fNeutronsCapturedPosZ.clear(); fNeutronsCapturedIsCaptureVolume.clear(); fNeutronsCapturedProductionE.clear(); + fNeutronsCapturedEDepByNeutron.clear(); + fNeutronsCapturedEDepByNeutronAndChildren.clear(); + fNeutronsCapturedEDepByNeutronInVeto.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVeto.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.clear(); fGammasNeutronCaptureNumber = 0; fGammasNeutronCapturePosX.clear(); @@ -278,10 +287,10 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv for (int j = 0; j < hits->GetNumberOfHits(); j++) { string process_name = (string)track->GetProcessName(hits->GetProcess(j)); if (process_name == "nCapture") { - cout << "Neutron capture!!!!!! " << particle_name << "trackId " << track->GetTrackID() - << " hit " << j << endl; + // << "Neutron capture!!!!!! " << particle_name << "trackId " << track->GetTrackID() + // << " hit " << j << endl; // track->PrintTrack(); - hits->PrintHits(j + 1); + // hits->PrintHits(j + 1); neutronsCaptured.insert(track->GetTrackID()); @@ -300,6 +309,55 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv } fNeutronsCapturedIsCaptureVolume.push_back(isCaptureVolume); fNeutronsCapturedProductionE.push_back(track->GetKineticEnergy()); + + // get energy deposited by neutron that undergoes capture and children + double neutronsCapturedEDepByNeutron = 0; + double neutronsCapturedEDepByNeutronAndChildren = 0; + double neutronsCapturedEDepByNeutronInVeto = 0; + double neutronsCapturedEDepByNeutronAndChildrenInVeto = 0; + + std::set parents = {track->GetTrackID()}; + std::map energy_in_veto; + for (int child = 0; child < fOutputG4Event->GetNumberOfTracks(); child++) { + auto track_child = fOutputG4Event->GetTrack(child); + if ((parents.count(track_child->GetParentID()) > 0) || + parents.count(track_child->GetTrackID()) > 0) { + // track or parent is in list of tracks, we add to list and add energy + parents.insert(track_child->GetTrackID()); + neutronsCapturedEDepByNeutronAndChildren += track_child->GetEnergy(); + if (track_child->GetTrackID() == track->GetTrackID()) { + neutronsCapturedEDepByNeutron += track_child->GetEnergy(); + } + for (const auto& vetoId : fVetoVolumeIds) { + neutronsCapturedEDepByNeutronAndChildrenInVeto += + track_child->GetEnergyInVolume(vetoId); + energy_in_veto[vetoId] += track_child->GetEnergyInVolume(vetoId); + if (track_child->GetTrackID() == track->GetTrackID()) { + neutronsCapturedEDepByNeutronInVeto += + track_child->GetEnergyInVolume(vetoId); + } + } + } + } + + fNeutronsCapturedEDepByNeutron.push_back(neutronsCapturedEDepByNeutron); + fNeutronsCapturedEDepByNeutronAndChildren.push_back( + neutronsCapturedEDepByNeutronAndChildren); + fNeutronsCapturedEDepByNeutronInVeto.push_back(neutronsCapturedEDepByNeutronInVeto); + fNeutronsCapturedEDepByNeutronAndChildrenInVeto.push_back( + neutronsCapturedEDepByNeutronAndChildrenInVeto); + + // get max and min energy in each veto (to compare with energy in ALL vetoes) + double energyMaxVeto = 0; + double energyMinVeto = -1; + for (const auto& pair : energy_in_veto) { + auto E = pair.second; + if (E > energyMaxVeto) energyMaxVeto = E; + if (E < energyMaxVeto || energyMinVeto == -1) energyMinVeto = E; + } + + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.push_back(energyMaxVeto); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.push_back(energyMinVeto); } } } @@ -311,7 +369,15 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv SetObservableValue("neutronsCapturedPosZ", fNeutronsCapturedPosZ); SetObservableValue("neutronsCapturedIsCaptureVolume", fNeutronsCapturedIsCaptureVolume); SetObservableValue("neutronsCapturedProductionE", fNeutronsCapturedProductionE); - + SetObservableValue("neutronsCapturedEDepByNeutron", fNeutronsCapturedEDepByNeutron); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildren", fNeutronsCapturedEDepByNeutronAndChildren); + SetObservableValue("neutronsCapturedEDepByNeutronInVeto", fNeutronsCapturedEDepByNeutronInVeto); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVeto", + fNeutronsCapturedEDepByNeutronAndChildrenInVeto); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMax", + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMin", + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin); for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { auto track = fOutputG4Event->GetTrack(i); string particle_name = (string)track->GetParticleName(); @@ -337,9 +403,9 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv fGammasNeutronCaptureIsCaptureVolume.push_back(isCaptureVolume); fGammasNeutronCaptureProductionE.push_back(track->GetKineticEnergy()); - cout << "gamma capture" << endl; + // cout << "gamma capture" << endl; - hits->PrintHits(1); + // hits->PrintHits(1); } } } @@ -404,15 +470,15 @@ TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEv } } - SetObservableValue("fSecondaryNeutronsShieldingNumber", fSecondaryNeutronsShieldingNumber); - SetObservableValue("fSecondaryNeutronsShieldingExitPosX", fSecondaryNeutronsShieldingExitPosX); - SetObservableValue("fSecondaryNeutronsShieldingExitPosY", fSecondaryNeutronsShieldingExitPosY); - SetObservableValue("fSecondaryNeutronsShieldingExitPosZ", fSecondaryNeutronsShieldingExitPosZ); - SetObservableValue("fSecondaryNeutronsShieldingIsCaptured", fSecondaryNeutronsShieldingIsCaptured); - SetObservableValue("fSecondaryNeutronsShieldingIsCapturedInCaptureVolume", + SetObservableValue("secondaryNeutronsShieldingNumber", fSecondaryNeutronsShieldingNumber); + SetObservableValue("secondaryNeutronsShieldingExitPosX", fSecondaryNeutronsShieldingExitPosX); + SetObservableValue("secondaryNeutronsShieldingExitPosY", fSecondaryNeutronsShieldingExitPosY); + SetObservableValue("secondaryNeutronsShieldingExitPosZ", fSecondaryNeutronsShieldingExitPosZ); + SetObservableValue("secondaryNeutronsShieldingIsCaptured", fSecondaryNeutronsShieldingIsCaptured); + SetObservableValue("secondaryNeutronsShieldingIsCapturedInCaptureVolume", fSecondaryNeutronsShieldingIsCapturedInCaptureVolume); - SetObservableValue("fSecondaryNeutronsShieldingProductionE", fSecondaryNeutronsShieldingProductionE); - SetObservableValue("fSecondaryNeutronsShieldingExitE", fSecondaryNeutronsShieldingExitE); + SetObservableValue("secondaryNeutronsShieldingProductionE", fSecondaryNeutronsShieldingProductionE); + SetObservableValue("secondaryNeutronsShieldingExitE", fSecondaryNeutronsShieldingExitE); return fOutputG4Event; } From a598009da2dd21a348d3895dab1a1057b6640395 Mon Sep 17 00:00:00 2001 From: lobis Date: Mon, 1 Mar 2021 13:37:45 +0100 Subject: [PATCH 4/4] updated docs --- ...eant4CosmicNeutronTaggingAnalysisProcess.h | 362 +++--- ...nt4CosmicNeutronTaggingAnalysisProcess.cxx | 1117 +++++++++-------- src/TRestGeant4VetoAnalysisProcess.cxx | 28 +- 3 files changed, 785 insertions(+), 722 deletions(-) diff --git a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h index 9791f0f2..e1bfc47d 100644 --- a/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h +++ b/inc/TRestGeant4CosmicNeutronTaggingAnalysisProcess.h @@ -1,181 +1,181 @@ -/************************************************************************* - * This file is part of the REST software framework. * - * * - * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * - * For more information see http://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 http://www.gnu.org/licenses/. * - * For the list of contributors see $REST_PATH/CREDITS. * - *************************************************************************/ - -#ifndef RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess -#define RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess - -#include -#include - -#include "TRestEventProcess.h" - -class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess { - private: - /// A pointer to the specific TRestGeant4Event input - TRestGeant4Event* fInputG4Event; //! - /// A pointer to the specific TRestGeant4Event output - TRestGeant4Event* fOutputG4Event; //! - /// A pointer to the simulation metadata information accessible to TRestRun - TRestGeant4Metadata* fG4Metadata; //! - - std::vector fVetoVolumeIds; //! - std::vector fCaptureVolumeIds; //! - std::vector fShieldingVolumeIds; //! - - string fVetoKeyword = ""; //! - string fCaptureKeyword = ""; //! - string fShieldingKeyword = ""; //! - std::vector fVetoGroupKeywords; //! - std::map> fVetoGroupVolumeNames; //! - std::vector fQuenchingFactors; //! - - // neutrons that undergo neutron capture - Int_t fNeutronsCapturedNumber; //! - std::vector fNeutronsCapturedPosX; //! - std::vector fNeutronsCapturedPosY; //! - std::vector fNeutronsCapturedPosZ; //! - std::vector fNeutronsCapturedIsCaptureVolume; //! - std::vector fNeutronsCapturedProductionE; //! - std::vector fNeutronsCapturedEDepByNeutron; //! - std::vector fNeutronsCapturedEDepByNeutronAndChildren; //! - std::vector fNeutronsCapturedEDepByNeutronInVeto; //! - std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVeto; //! - std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax; //! - std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin; //! - - // gammas that are produced from neutron capture - Int_t fGammasNeutronCaptureNumber; //! - std::vector fGammasNeutronCapturePosX; //! - std::vector fGammasNeutronCapturePosY; //! - std::vector fGammasNeutronCapturePosZ; //! - std::vector fGammasNeutronCaptureIsCaptureVolume; //! - std::vector fGammasNeutronCaptureProductionE; //! - - // secondary neutrons that exit the lead shielding - Int_t fSecondaryNeutronsShieldingNumber; //! - std::vector fSecondaryNeutronsShieldingExitPosX; //! - std::vector fSecondaryNeutronsShieldingExitPosY; //! - std::vector fSecondaryNeutronsShieldingExitPosZ; //! - std::vector fSecondaryNeutronsShieldingIsCaptured; //! - std::vector fSecondaryNeutronsShieldingIsCapturedInCaptureVolume; //! - std::vector fSecondaryNeutronsShieldingProductionE; //! - std::vector fSecondaryNeutronsShieldingExitE; //! - - void InitFromConfigFile(); - void Initialize(); - void LoadDefaultConfig(); - void Reset(); - // clean string (https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring) - inline std::string& rtrim(std::string& s, const char* t = " \t\n\r\f\v") { - s.erase(s.find_last_not_of(t) + 1); - return s; - } - // trim from beginning of string (left) - inline std::string& ltrim(std::string& s, const char* t = " \t\n\r\f\v") { - s.erase(0, s.find_first_not_of(t)); - return s; - } - // trim from both ends of string (right then left) - inline std::string& trim(std::string& s, const char* t = " \t\n\r\f\v") { return ltrim(rtrim(s, t), t); } - - // final clean string: trim and UPPER - inline std::string& clean_string(std::string& s) { - s = trim(s); - std::transform(s.begin(), s.end(), s.begin(), ::tolower); - return s; - } - - protected: - // add here the members of your event process - - public: - any GetInputEvent() { return fInputG4Event; } - any GetOutputEvent() { return fOutputG4Event; } - - void InitProcess(); - TRestEvent* ProcessEvent(TRestEvent* eventInput); - void EndProcess(); - - void LoadConfig(std::string cfgFilename, std::string name = ""); - - /// It prints out the process parameters stored in the metadata structure - void PrintMetadata() { - BeginPrintProcess(); - - debug << "VETO KEYWORD: " << fVetoKeyword << endl; - debug << endl; - - debug << "VETO GROUP KEYWORDS:" << endl; - for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { - debug << "\t" << fVetoGroupKeywords[i] << endl; - } - debug << endl; - - debug << "Found " << fVetoVolumeIds.size() << " veto volumes:" << endl; - for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { - debug << "\t" << fG4Metadata->GetActiveVolumeName(fVetoVolumeIds[i]) << endl; - } - debug << endl; - - // capture volumes - - debug << "CAPTURE KEYWORD: " << fCaptureKeyword << endl; - debug << endl; - - debug << "Found " << fCaptureVolumeIds.size() << " Capture volumes:" << endl; - for (unsigned int i = 0; i < fCaptureVolumeIds.size(); i++) { - debug << "\t" << fG4Metadata->GetActiveVolumeName(fCaptureVolumeIds[i]) << endl; - } - debug << endl; - - // shielding volume/s - - debug << "SHIELDING KEYWORD: " << fShieldingKeyword << endl; - debug << endl; - - debug << "Found " << fShieldingVolumeIds.size() << " Shielding volumes:" << endl; - for (unsigned int i = 0; i < fShieldingVolumeIds.size(); i++) { - debug << "\t" << fG4Metadata->GetActiveVolumeName(fShieldingVolumeIds[i]) << endl; - } - debug << endl; - - debug << "QUENCHING FACTORS (" << fQuenchingFactors.size() << " Total)" << endl; - for (unsigned int i = 0; i < fQuenchingFactors.size(); i++) { - debug << "\t" << fQuenchingFactors[i] << endl; - } - debug << endl; - - EndPrintProcess(); - } - - /// Returns a new instance of this class - TRestEventProcess* Maker() { return new TRestGeant4CosmicNeutronTaggingAnalysisProcess; } - /// Returns the name of this process - TString GetProcessName() { return (TString) "geant4VetoAnalysis"; } - - TRestGeant4CosmicNeutronTaggingAnalysisProcess(); - TRestGeant4CosmicNeutronTaggingAnalysisProcess(char* cfgFileName); - ~TRestGeant4CosmicNeutronTaggingAnalysisProcess(); - - ClassDef(TRestGeant4CosmicNeutronTaggingAnalysisProcess, 1); -}; -#endif // RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +#ifndef RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess +#define RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess + +#include +#include + +#include "TRestEventProcess.h" + +class TRestGeant4CosmicNeutronTaggingAnalysisProcess : public TRestEventProcess { + private: + /// A pointer to the specific TRestGeant4Event input + TRestGeant4Event* fInputG4Event; //! + /// A pointer to the specific TRestGeant4Event output + TRestGeant4Event* fOutputG4Event; //! + /// A pointer to the simulation metadata information accessible to TRestRun + TRestGeant4Metadata* fG4Metadata; //! + + std::vector fVetoVolumeIds; //! + std::vector fCaptureVolumeIds; //! + std::vector fShieldingVolumeIds; //! + + string fVetoKeyword = ""; //! + string fCaptureKeyword = ""; //! + string fShieldingKeyword = ""; //! + std::vector fVetoGroupKeywords; //! + std::map> fVetoGroupVolumeNames; //! + std::vector fQuenchingFactors; //! + + // neutrons that undergo neutron capture + Int_t fNeutronsCapturedNumber; //! + std::vector fNeutronsCapturedPosX; //! + std::vector fNeutronsCapturedPosY; //! + std::vector fNeutronsCapturedPosZ; //! + std::vector fNeutronsCapturedIsCaptureVolume; //! + std::vector fNeutronsCapturedProductionE; //! + std::vector fNeutronsCapturedEDepByNeutron; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildren; //! + std::vector fNeutronsCapturedEDepByNeutronInVeto; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVeto; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax; //! + std::vector fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin; //! + + // gammas that are produced from neutron capture + Int_t fGammasNeutronCaptureNumber; //! + std::vector fGammasNeutronCapturePosX; //! + std::vector fGammasNeutronCapturePosY; //! + std::vector fGammasNeutronCapturePosZ; //! + std::vector fGammasNeutronCaptureIsCaptureVolume; //! + std::vector fGammasNeutronCaptureProductionE; //! + + // secondary neutrons that exit the lead shielding + Int_t fSecondaryNeutronsShieldingNumber; //! + std::vector fSecondaryNeutronsShieldingExitPosX; //! + std::vector fSecondaryNeutronsShieldingExitPosY; //! + std::vector fSecondaryNeutronsShieldingExitPosZ; //! + std::vector fSecondaryNeutronsShieldingIsCaptured; //! + std::vector fSecondaryNeutronsShieldingIsCapturedInCaptureVolume; //! + std::vector fSecondaryNeutronsShieldingProductionE; //! + std::vector fSecondaryNeutronsShieldingExitE; //! + + void InitFromConfigFile(); + void Initialize(); + void LoadDefaultConfig(); + void Reset(); + // clean string (https://stackoverflow.com/questions/216823/whats-the-best-way-to-trim-stdstring) + inline std::string& rtrim(std::string& s, const char* t = " \t\n\r\f\v") { + s.erase(s.find_last_not_of(t) + 1); + return s; + } + // trim from beginning of string (left) + inline std::string& ltrim(std::string& s, const char* t = " \t\n\r\f\v") { + s.erase(0, s.find_first_not_of(t)); + return s; + } + // trim from both ends of string (right then left) + inline std::string& trim(std::string& s, const char* t = " \t\n\r\f\v") { return ltrim(rtrim(s, t), t); } + + // final clean string: trim and UPPER + inline std::string& clean_string(std::string& s) { + s = trim(s); + std::transform(s.begin(), s.end(), s.begin(), ::tolower); + return s; + } + + protected: + // add here the members of your event process + + public: + any GetInputEvent() { return fInputG4Event; } + any GetOutputEvent() { return fOutputG4Event; } + + void InitProcess(); + TRestEvent* ProcessEvent(TRestEvent* eventInput); + void EndProcess(); + + void LoadConfig(std::string cfgFilename, std::string name = ""); + + /// It prints out the process parameters stored in the metadata structure + void PrintMetadata() { + BeginPrintProcess(); + + debug << "VETO KEYWORD: " << fVetoKeyword << endl; + debug << endl; + + debug << "VETO GROUP KEYWORDS:" << endl; + for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { + debug << "\t" << fVetoGroupKeywords[i] << endl; + } + debug << endl; + + debug << "Found " << fVetoVolumeIds.size() << " veto volumes:" << endl; + for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fVetoVolumeIds[i]) << endl; + } + debug << endl; + + // capture volumes + + debug << "CAPTURE KEYWORD: " << fCaptureKeyword << endl; + debug << endl; + + debug << "Found " << fCaptureVolumeIds.size() << " Capture volumes:" << endl; + for (unsigned int i = 0; i < fCaptureVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fCaptureVolumeIds[i]) << endl; + } + debug << endl; + + // shielding volume/s + + debug << "SHIELDING KEYWORD: " << fShieldingKeyword << endl; + debug << endl; + + debug << "Found " << fShieldingVolumeIds.size() << " Shielding volumes:" << endl; + for (unsigned int i = 0; i < fShieldingVolumeIds.size(); i++) { + debug << "\t" << fG4Metadata->GetActiveVolumeName(fShieldingVolumeIds[i]) << endl; + } + debug << endl; + + debug << "QUENCHING FACTORS (" << fQuenchingFactors.size() << " Total)" << endl; + for (unsigned int i = 0; i < fQuenchingFactors.size(); i++) { + debug << "\t" << fQuenchingFactors[i] << endl; + } + debug << endl; + + EndPrintProcess(); + } + + /// Returns a new instance of this class + TRestEventProcess* Maker() { return new TRestGeant4CosmicNeutronTaggingAnalysisProcess; } + /// Returns the name of this process + TString GetProcessName() { return (TString) "geant4VetoAnalysis"; } + + TRestGeant4CosmicNeutronTaggingAnalysisProcess(); + TRestGeant4CosmicNeutronTaggingAnalysisProcess(char* cfgFileName); + ~TRestGeant4CosmicNeutronTaggingAnalysisProcess(); + + ClassDef(TRestGeant4CosmicNeutronTaggingAnalysisProcess, 1); +}; +#endif // RestCore_TRestGeant4CosmicNeutronTaggingAnalysisProcess diff --git a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx index 01c7a9ca..3470727c 100644 --- a/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx +++ b/src/TRestGeant4CosmicNeutronTaggingAnalysisProcess.cxx @@ -1,540 +1,577 @@ -/************************************************************************* - * This file is part of the REST software framework. * - * * - * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * - * For more information see http://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 http://www.gnu.org/licenses/. * - * For the list of contributors see $REST_PATH/CREDITS. * - *************************************************************************/ - -////////////////////////////////////////////////////////////////////////// -/// TODO: fill documentation -/// -///-------------------------------------------------------------------------- -/// -/// RESTsoft - Software for Rare Event Searches with TPCs -/// -/// History of developments: -/// -/// 2021-February: Implementation. -/// -/// \class TRestGeant4CosmicNeutronTaggingAnalysisProcess -/// \author Luis Obis -/// -///
-/// - -#include "TRestGeant4CosmicNeutronTaggingAnalysisProcess.h" -using namespace std; - -ClassImp(TRestGeant4CosmicNeutronTaggingAnalysisProcess); - -TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess() { - Initialize(); -} - -TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess( - char* cfgFileName) { - Initialize(); - if (LoadConfigFromFile(cfgFileName)) LoadDefaultConfig(); -} - -/////////////////////////////////////////////// -/// \brief Default destructor -/// -TRestGeant4CosmicNeutronTaggingAnalysisProcess::~TRestGeant4CosmicNeutronTaggingAnalysisProcess() { - delete fOutputG4Event; -} - -/////////////////////////////////////////////// -/// \brief Function to load the default config in absence of RML input -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); } - -/////////////////////////////////////////////// -/// \brief Function to initialize input/output event members and define the -/// section name -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Initialize() { - fG4Metadata = NULL; - - SetSectionName(this->ClassName()); - SetLibraryVersion(LIBRARY_VERSION); - - fInputG4Event = NULL; - fOutputG4Event = new TRestGeant4Event(); -} - -/////////////////////////////////////////////// -/// \brief Function to load the configuration from an external configuration -/// file. -/// -/// If no configuration path is defined in TRestMetadata::SetConfigFilePath -/// the path to the config file must be specified using full path, absolute or -/// relative. -/// -/// \param cfgFileName A const char* giving the path to an RML file. -/// \param name The name of the specific metadata. It will be used to find the -/// correspondig TRestGeant4CosmicNeutronTaggingAnalysisProcess section inside the RML. -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadConfig(std::string cfgFilename, std::string name) { - if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); -} - -/////////////////////////////////////////////// -/// \brief Process initialization. -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { - fG4Metadata = GetMetadata(); - - // CAREFUL THIS METHOD IS CALLED TWICE! - // we need to reset these variables to zero - Reset(); - // get "veto" volumes - if (fVetoVolumeIds.empty()) { - for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); - volume_name = clean_string(volume_name); - if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { - fVetoVolumeIds.push_back(i); - } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { - fCaptureVolumeIds.push_back(i); - } else if (volume_name.find(clean_string(fShieldingKeyword)) != string::npos) { - fShieldingVolumeIds.push_back(i); - } - } - - // veto groups (fill fVetoGroupVolumeNames) - for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { - string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); - fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; - for (int& id : fVetoVolumeIds) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); - volume_name = clean_string(volume_name); - if (volume_name.find(veto_group_keyword) != string::npos) { - fVetoGroupVolumeNames[veto_group_keyword].push_back( - (string)fG4Metadata->GetActiveVolumeName(id)); - } - } - } - } - - PrintMetadata(); -} - -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Reset() { - /* - fVetoVolumeIds.clear(); - fVetoGroupVolumeNames.clear(); - fCaptureVolumeIds.clear(); - */ - fNeutronsCapturedNumber = 0; - fNeutronsCapturedPosX.clear(); - fNeutronsCapturedPosY.clear(); - fNeutronsCapturedPosZ.clear(); - fNeutronsCapturedIsCaptureVolume.clear(); - fNeutronsCapturedProductionE.clear(); - fNeutronsCapturedEDepByNeutron.clear(); - fNeutronsCapturedEDepByNeutronAndChildren.clear(); - fNeutronsCapturedEDepByNeutronInVeto.clear(); - fNeutronsCapturedEDepByNeutronAndChildrenInVeto.clear(); - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.clear(); - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.clear(); - - fGammasNeutronCaptureNumber = 0; - fGammasNeutronCapturePosX.clear(); - fGammasNeutronCapturePosY.clear(); - fGammasNeutronCapturePosZ.clear(); - fGammasNeutronCaptureIsCaptureVolume.clear(); - fGammasNeutronCaptureProductionE.clear(); - - fSecondaryNeutronsShieldingNumber = 0; - fSecondaryNeutronsShieldingExitPosX.clear(); - fSecondaryNeutronsShieldingExitPosY.clear(); - fSecondaryNeutronsShieldingExitPosZ.clear(); - fSecondaryNeutronsShieldingIsCaptured.clear(); - fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.clear(); - fSecondaryNeutronsShieldingProductionE.clear(); - fSecondaryNeutronsShieldingExitE.clear(); -} -/////////////////////////////////////////////// -/// \brief The main processing event function -/// -TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEvent* evInput) { - fInputG4Event = (TRestGeant4Event*)evInput; - *fOutputG4Event = *((TRestGeant4Event*)evInput); - - Reset(); - std::map volume_energy_map; - - for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { - int id = fVetoVolumeIds[i]; - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); - - Double_t energy = fOutputG4Event->GetEnergyDepositedInVolume(id); - volume_energy_map[volume_name] = energy; - } - - Double_t energy_veto_max = 0; - for (const auto& pair : volume_energy_map) { - Double_t veto_energy = pair.second; - SetObservableValue(pair.first + "VolumeEDep", veto_energy); - if (veto_energy > energy_veto_max) { - energy_veto_max = veto_energy; - }; - } - SetObservableValue("vetoAllEVetoMax", energy_veto_max); - - // veto groups - for (const auto& pair : fVetoGroupVolumeNames) { - Double_t energy_veto_max_group = 0; - for (unsigned int i = 0; i < pair.second.size(); i++) { - string volume_name = pair.second[i]; - Double_t veto_energy = volume_energy_map[volume_name]; - if (veto_energy > energy_veto_max_group) { - energy_veto_max_group = veto_energy; - }; - } - // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) - string group_name; - for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { - if (it == pair.first.cbegin()) { - group_name += std::toupper(*it); - } else { - group_name += std::tolower(*it); - } - } - SetObservableValue("vetoGroup" + group_name + "EVetoMax", energy_veto_max_group); - } - - // - - for (const auto& quenching_factor : fQuenchingFactors) { - string quenching_factor_string = std::to_string(quenching_factor); - // replace "." in string by "_" because its gives very strange problems - quenching_factor_string = - quenching_factor_string.replace(quenching_factor_string.find("."), sizeof(".") - 1, "_"); - volume_energy_map.clear(); - for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { - auto track = fOutputG4Event->GetTrack(i); - string particle_name = (string)track->GetParticleName(); - for (const auto& id : fVetoVolumeIds) { - string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); - - if (particle_name == "e-" || particle_name == "e+" || particle_name == "gamma") { - // no quenching factor - volume_energy_map[volume_name] += track->GetEnergyInVolume(id); - } else { - // apply quenching factor - volume_energy_map[volume_name] += quenching_factor * track->GetEnergyInVolume(id); - } - } - } - - Double_t energy_veto_max = 0; - for (const auto& pair : volume_energy_map) { - Double_t veto_energy = pair.second; - SetObservableValue(pair.first + "VolumeEDep" + "Qf" + quenching_factor_string, veto_energy); - if (veto_energy > energy_veto_max) { - energy_veto_max = veto_energy; - }; - } - SetObservableValue(string("vetoAllEVetoMax") + "Qf" + quenching_factor_string, energy_veto_max); - - // veto groups - for (const auto& pair : fVetoGroupVolumeNames) { - Double_t energy_veto_max_group = 0; - for (unsigned int i = 0; i < pair.second.size(); i++) { - string volume_name = pair.second[i]; - Double_t veto_energy = volume_energy_map[volume_name]; - if (veto_energy > energy_veto_max_group) { - energy_veto_max_group = veto_energy; - }; - } - // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) - string group_name; - for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { - if (it == pair.first.cbegin()) { - group_name += std::toupper(*it); - } else { - group_name += std::tolower(*it); - } - } - SetObservableValue("vetoGroup" + group_name + "EVetoMax" + "Qf" + quenching_factor_string, - energy_veto_max_group); - } - } - - std::set neutronsCaptured = {}; - for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { - auto track = fOutputG4Event->GetTrack(i); - string particle_name = (string)track->GetParticleName(); - if (particle_name == "neutron") { - auto hits = track->GetHits(); - for (int j = 0; j < hits->GetNumberOfHits(); j++) { - string process_name = (string)track->GetProcessName(hits->GetProcess(j)); - if (process_name == "nCapture") { - // << "Neutron capture!!!!!! " << particle_name << "trackId " << track->GetTrackID() - // << " hit " << j << endl; - // track->PrintTrack(); - // hits->PrintHits(j + 1); - - neutronsCaptured.insert(track->GetTrackID()); - - fNeutronsCapturedNumber += 1; - fNeutronsCapturedPosX.push_back(hits->GetX(j)); - fNeutronsCapturedPosY.push_back(hits->GetY(j)); - fNeutronsCapturedPosZ.push_back(hits->GetZ(j)); - - Int_t volumeId = hits->GetVolumeId(j); - Int_t isCaptureVolume = 0; - for (const auto& id : fCaptureVolumeIds) { - if (volumeId == id) { - isCaptureVolume = 1; - continue; - } - } - fNeutronsCapturedIsCaptureVolume.push_back(isCaptureVolume); - fNeutronsCapturedProductionE.push_back(track->GetKineticEnergy()); - - // get energy deposited by neutron that undergoes capture and children - double neutronsCapturedEDepByNeutron = 0; - double neutronsCapturedEDepByNeutronAndChildren = 0; - double neutronsCapturedEDepByNeutronInVeto = 0; - double neutronsCapturedEDepByNeutronAndChildrenInVeto = 0; - - std::set parents = {track->GetTrackID()}; - std::map energy_in_veto; - for (int child = 0; child < fOutputG4Event->GetNumberOfTracks(); child++) { - auto track_child = fOutputG4Event->GetTrack(child); - if ((parents.count(track_child->GetParentID()) > 0) || - parents.count(track_child->GetTrackID()) > 0) { - // track or parent is in list of tracks, we add to list and add energy - parents.insert(track_child->GetTrackID()); - neutronsCapturedEDepByNeutronAndChildren += track_child->GetEnergy(); - if (track_child->GetTrackID() == track->GetTrackID()) { - neutronsCapturedEDepByNeutron += track_child->GetEnergy(); - } - for (const auto& vetoId : fVetoVolumeIds) { - neutronsCapturedEDepByNeutronAndChildrenInVeto += - track_child->GetEnergyInVolume(vetoId); - energy_in_veto[vetoId] += track_child->GetEnergyInVolume(vetoId); - if (track_child->GetTrackID() == track->GetTrackID()) { - neutronsCapturedEDepByNeutronInVeto += - track_child->GetEnergyInVolume(vetoId); - } - } - } - } - - fNeutronsCapturedEDepByNeutron.push_back(neutronsCapturedEDepByNeutron); - fNeutronsCapturedEDepByNeutronAndChildren.push_back( - neutronsCapturedEDepByNeutronAndChildren); - fNeutronsCapturedEDepByNeutronInVeto.push_back(neutronsCapturedEDepByNeutronInVeto); - fNeutronsCapturedEDepByNeutronAndChildrenInVeto.push_back( - neutronsCapturedEDepByNeutronAndChildrenInVeto); - - // get max and min energy in each veto (to compare with energy in ALL vetoes) - double energyMaxVeto = 0; - double energyMinVeto = -1; - for (const auto& pair : energy_in_veto) { - auto E = pair.second; - if (E > energyMaxVeto) energyMaxVeto = E; - if (E < energyMaxVeto || energyMinVeto == -1) energyMinVeto = E; - } - - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.push_back(energyMaxVeto); - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.push_back(energyMinVeto); - } - } - } - } - - SetObservableValue("neutronsCapturedNumber", fNeutronsCapturedNumber); - SetObservableValue("neutronsCapturedPosX", fNeutronsCapturedPosX); - SetObservableValue("neutronsCapturedPosY", fNeutronsCapturedPosY); - SetObservableValue("neutronsCapturedPosZ", fNeutronsCapturedPosZ); - SetObservableValue("neutronsCapturedIsCaptureVolume", fNeutronsCapturedIsCaptureVolume); - SetObservableValue("neutronsCapturedProductionE", fNeutronsCapturedProductionE); - SetObservableValue("neutronsCapturedEDepByNeutron", fNeutronsCapturedEDepByNeutron); - SetObservableValue("neutronsCapturedEDepByNeutronAndChildren", fNeutronsCapturedEDepByNeutronAndChildren); - SetObservableValue("neutronsCapturedEDepByNeutronInVeto", fNeutronsCapturedEDepByNeutronInVeto); - SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVeto", - fNeutronsCapturedEDepByNeutronAndChildrenInVeto); - SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMax", - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax); - SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMin", - fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin); - for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { - auto track = fOutputG4Event->GetTrack(i); - string particle_name = (string)track->GetParticleName(); - if (particle_name == "gamma") { - // check if gamma is child of captured neutron - Int_t parent = track->GetParentID(); - if (neutronsCaptured.count(parent) > 0) { - auto hits = track->GetHits(); - - fGammasNeutronCaptureNumber += 1; - fGammasNeutronCapturePosX.push_back(hits->GetX(0)); - fGammasNeutronCapturePosY.push_back(hits->GetY(0)); - fGammasNeutronCapturePosZ.push_back(hits->GetZ(0)); - - Int_t volumeId = hits->GetVolumeId(0); - Int_t isCaptureVolume = 0; - for (const auto& id : fCaptureVolumeIds) { - if (volumeId == id) { - isCaptureVolume = 1; - continue; - } - } - fGammasNeutronCaptureIsCaptureVolume.push_back(isCaptureVolume); - fGammasNeutronCaptureProductionE.push_back(track->GetKineticEnergy()); - - // cout << "gamma capture" << endl; - - // hits->PrintHits(1); - } - } - } - - SetObservableValue("gammasNeutronCaptureNumber", fGammasNeutronCaptureNumber); - SetObservableValue("gammasNeutronCapturePosX", fGammasNeutronCapturePosX); - SetObservableValue("gammasNeutronCapturePosY", fGammasNeutronCapturePosY); - SetObservableValue("gammasNeutronCapturePosZ", fGammasNeutronCapturePosZ); - SetObservableValue("gammasNeutronCaptureIsCaptureVolume", fGammasNeutronCaptureIsCaptureVolume); - SetObservableValue("gammasNeutronCaptureProductionE", fGammasNeutronCaptureProductionE); - - std::set secondaryNeutrons = {}; // avoid counting twice - for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { - auto track = fOutputG4Event->GetTrack(i); - string particle_name = (string)track->GetParticleName(); - if (particle_name == "neutron" && track->GetParentID() != 0) { // not consider primary - // check if neutron exits shielding - auto hits = track->GetHits(); - for (int j = 0; j < hits->GetNumberOfHits(); j++) { - string process_name = (string)track->GetProcessName(hits->GetProcess(j)); - if (process_name == "Transportation") { - for (const auto& id : fShieldingVolumeIds) { - if (hits->GetVolumeId(j) == id) { - // transportation and shielding == exits shielding - if (secondaryNeutrons.count(track->GetTrackID()) == 0) { - // first time adding this secondary neutron - secondaryNeutrons.insert(track->GetTrackID()); - } else { - continue; - } - fSecondaryNeutronsShieldingNumber += 1; - fSecondaryNeutronsShieldingExitPosX.push_back(hits->GetX(j)); - fSecondaryNeutronsShieldingExitPosY.push_back(hits->GetY(j)); - fSecondaryNeutronsShieldingExitPosZ.push_back(hits->GetZ(j)); - - Int_t volumeId = hits->GetVolumeId(j); - Int_t isCaptureVolume = 0; - for (const auto& id : fCaptureVolumeIds) { - if (volumeId == id) { - isCaptureVolume = 1; - continue; - } - } - Int_t isCaptured = 0; - if (neutronsCaptured.count(track->GetTrackID()) > 0) { - isCaptured = 1; - } - fSecondaryNeutronsShieldingIsCaptured.push_back(isCaptured); - if (isCaptured) - fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back( - isCaptureVolume); - else { - fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back(0); - } - - fSecondaryNeutronsShieldingProductionE.push_back(track->GetKineticEnergy()); - fSecondaryNeutronsShieldingExitE.push_back(hits->GetKineticEnergy(j)); - } - } - } - } - } - } - - SetObservableValue("secondaryNeutronsShieldingNumber", fSecondaryNeutronsShieldingNumber); - SetObservableValue("secondaryNeutronsShieldingExitPosX", fSecondaryNeutronsShieldingExitPosX); - SetObservableValue("secondaryNeutronsShieldingExitPosY", fSecondaryNeutronsShieldingExitPosY); - SetObservableValue("secondaryNeutronsShieldingExitPosZ", fSecondaryNeutronsShieldingExitPosZ); - SetObservableValue("secondaryNeutronsShieldingIsCaptured", fSecondaryNeutronsShieldingIsCaptured); - SetObservableValue("secondaryNeutronsShieldingIsCapturedInCaptureVolume", - fSecondaryNeutronsShieldingIsCapturedInCaptureVolume); - SetObservableValue("secondaryNeutronsShieldingProductionE", fSecondaryNeutronsShieldingProductionE); - SetObservableValue("secondaryNeutronsShieldingExitE", fSecondaryNeutronsShieldingExitE); - - return fOutputG4Event; -} - -/////////////////////////////////////////////// -/// \brief Function to include required actions after all events have been -/// processed. -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::EndProcess() { - // Function to be executed once at the end of the process - // (after all events have been processed) - - // Start by calling the EndProcess function of the abstract class. - // Comment this if you don't want it. - // TRestEventProcess::EndProcess(); -} - -/////////////////////////////////////////////// -/// \brief Function to read input parameters from the RML -/// TRestGeant4CosmicNeutronTaggingAnalysisProcess metadata section -/// -void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitFromConfigFile() { - // word to identify active volume as veto (default = "veto" e.g. "vetoTop") - string veto_keyword = GetParameter("vetoKeyword", "veto"); - fVetoKeyword = clean_string(veto_keyword); - // comma separated tags: "top, bottom, ..." - string veto_group_keywords = GetParameter("vetoGroupKeywords", ""); - stringstream ss(veto_group_keywords); - while (ss.good()) { - string substr; - getline(ss, substr, ','); - fVetoGroupKeywords.push_back(clean_string(substr)); - } - - // word to identify active volume as capture sheet (cadmium, default = "sheet" e.g. - // "scintillatorSheetTop1of4") - string capture_keyword = GetParameter("captureKeyword", "sheet"); - fCaptureKeyword = clean_string(capture_keyword); - - // word to identify active volume as shielding - - string shielding_keyword = GetParameter("shieldingKeyword", "shielding"); - fShieldingKeyword = clean_string(shielding_keyword); - - // comma separated quenching factors: "0.15, 1.00, ..." - string quenching_factors = GetParameter("vetoQuenchingFactors", "-1"); - stringstream ss_qf(quenching_factors); - while (ss_qf.good()) { - string substr; - getline(ss_qf, substr, ','); - substr = clean_string(substr); - Float_t quenching_factor = (Float_t)std::atof(substr.c_str()); - if (quenching_factor > 1 || quenching_factor < 0) { - cout << "ERROR: quenching factor must be between 0 and 1" << endl; - continue; - } - fQuenchingFactors.push_back(quenching_factor); - } -} +/************************************************************************* + * This file is part of the REST software framework. * + * * + * Copyright (C) 2020 GIFNA/TREX (University of Zaragoza) * + * For more information see http://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 http://www.gnu.org/licenses/. * + * For the list of contributors see $REST_PATH/CREDITS. * + *************************************************************************/ + +////////////////////////////////////////////////////////////////////////// +/// The TRestGeant4CosmicNeutronTaggingAnalysisProcess generates observables based on veto volumes energy +/// depositions. It was first developed as a process for the IAXO experiment but can be used in any analysis. +/// It uses `keywords` to identify different relevant volumes (such as vetoes). The default veto keyword for +/// IAXO is `veto` and it will tag each volume containing the keyword as a veto volume, so avoid using the +/// keyword on volumes that do not act as vetoes. There are also keywords for shielding and capture volumes +/// (Cd layers). +/// +/// Parameters: +/// `vetoKeyword`: keyword to identify a volume of the geometry as a veto. The condition is that the keyword +/// is contained inside the name of the volume. Only volumes serving as veto (i.e. scintillators) should +/// contain this vetoKeyword in their names. +/// `captureKeyword`: keyword to identify a volume of the geometry as a capture volumes. These volumes +/// correspond to the volumes used to capture neutrons and produce easily detectable secondaries, such as Cd +/// layers. This parameter is optional and is useful to benchmark the effect of the capture volume and +/// material, for example, if a very low % of neutron captures occur on the `capture` volumes, they are not +/// very optimized. +/// `shieldingKeyword`: keyword to identify the shielding volume. It is used to study the secondaries coming +/// out of the shielding, as in IAXO most of the secondaries come from the shielding. If there are multiple +/// shielding volumes this may not work as expected. +/// `vetoGroupKeywords`: comma separated keywords used to identify different groups of vetoes. This is an +/// optional parameter that when specified will make the process return additional observables on a per group +/// basis. The most common use case is using group names to identify the location of the vetoes (top, bottom, +/// ...). The volumes detected as vetoes (via vetoKeyword) will also be assigned to a group if they contain a +/// keyword contained in this list. +/// `vetoQuenchingFactors`: comma separated values for the quenching factors used in the analysis. The +/// observables will be calculated for each of the quenching factors contained in this list. Values between 0 +/// and 1 only. This is useful in case the user doesn't know the exact value of the quenching factor. Also it +/// is useful to insert `0` or `1` to study the effects of electromagnetic processes only, or no quenching at +/// all. +/// +/// Example usage: +/// +/// +/// +/// +/// +/// +/// +///-------------------------------------------------------------------------- +/// +/// RESTsoft - Software for Rare Event Searches with TPCs +/// +/// History of developments: +/// +/// 2021-February: Implementation. +/// +/// \class TRestGeant4CosmicNeutronTaggingAnalysisProcess +/// \author Luis Obis +/// +///
+/// + +#include "TRestGeant4CosmicNeutronTaggingAnalysisProcess.h" +using namespace std; + +ClassImp(TRestGeant4CosmicNeutronTaggingAnalysisProcess); + +TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess() { + Initialize(); +} + +TRestGeant4CosmicNeutronTaggingAnalysisProcess::TRestGeant4CosmicNeutronTaggingAnalysisProcess( + char* cfgFileName) { + Initialize(); + if (LoadConfigFromFile(cfgFileName)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Default destructor +/// +TRestGeant4CosmicNeutronTaggingAnalysisProcess::~TRestGeant4CosmicNeutronTaggingAnalysisProcess() { + delete fOutputG4Event; +} + +/////////////////////////////////////////////// +/// \brief Function to load the default config in absence of RML input +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadDefaultConfig() { SetTitle("Default config"); } + +/////////////////////////////////////////////// +/// \brief Function to initialize input/output event members and define the +/// section name +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Initialize() { + fG4Metadata = NULL; + + SetSectionName(this->ClassName()); + SetLibraryVersion(LIBRARY_VERSION); + + fInputG4Event = NULL; + fOutputG4Event = new TRestGeant4Event(); +} + +/////////////////////////////////////////////// +/// \brief Function to load the configuration from an external configuration +/// file. +/// +/// If no configuration path is defined in TRestMetadata::SetConfigFilePath +/// the path to the config file must be specified using full path, absolute or +/// relative. +/// +/// \param cfgFileName A const char* giving the path to an RML file. +/// \param name The name of the specific metadata. It will be used to find the +/// correspondig TRestGeant4CosmicNeutronTaggingAnalysisProcess section inside the RML. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::LoadConfig(std::string cfgFilename, std::string name) { + if (LoadConfigFromFile(cfgFilename, name)) LoadDefaultConfig(); +} + +/////////////////////////////////////////////// +/// \brief Process initialization. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitProcess() { + fG4Metadata = GetMetadata(); + + // CAREFUL THIS METHOD IS CALLED TWICE! + // we need to reset these variables to zero + Reset(); + // get "veto" volumes + if (fVetoVolumeIds.empty()) { + for (unsigned int i = 0; i < fG4Metadata->GetNumberOfActiveVolumes(); i++) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(i); + volume_name = clean_string(volume_name); + if (volume_name.find(clean_string(fVetoKeyword)) != string::npos) { + fVetoVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fCaptureKeyword)) != string::npos) { + fCaptureVolumeIds.push_back(i); + } else if (volume_name.find(clean_string(fShieldingKeyword)) != string::npos) { + fShieldingVolumeIds.push_back(i); + } + } + + // veto groups (fill fVetoGroupVolumeNames) + for (unsigned int i = 0; i < fVetoGroupKeywords.size(); i++) { + string veto_group_keyword = clean_string(fVetoGroupKeywords[i]); + fVetoGroupVolumeNames[veto_group_keyword] = std::vector{}; + for (int& id : fVetoVolumeIds) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + volume_name = clean_string(volume_name); + if (volume_name.find(veto_group_keyword) != string::npos) { + fVetoGroupVolumeNames[veto_group_keyword].push_back( + (string)fG4Metadata->GetActiveVolumeName(id)); + } + } + } + } + + PrintMetadata(); +} + +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::Reset() { + /* + fVetoVolumeIds.clear(); + fVetoGroupVolumeNames.clear(); + fCaptureVolumeIds.clear(); + */ + fNeutronsCapturedNumber = 0; + fNeutronsCapturedPosX.clear(); + fNeutronsCapturedPosY.clear(); + fNeutronsCapturedPosZ.clear(); + fNeutronsCapturedIsCaptureVolume.clear(); + fNeutronsCapturedProductionE.clear(); + fNeutronsCapturedEDepByNeutron.clear(); + fNeutronsCapturedEDepByNeutronAndChildren.clear(); + fNeutronsCapturedEDepByNeutronInVeto.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVeto.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.clear(); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.clear(); + + fGammasNeutronCaptureNumber = 0; + fGammasNeutronCapturePosX.clear(); + fGammasNeutronCapturePosY.clear(); + fGammasNeutronCapturePosZ.clear(); + fGammasNeutronCaptureIsCaptureVolume.clear(); + fGammasNeutronCaptureProductionE.clear(); + + fSecondaryNeutronsShieldingNumber = 0; + fSecondaryNeutronsShieldingExitPosX.clear(); + fSecondaryNeutronsShieldingExitPosY.clear(); + fSecondaryNeutronsShieldingExitPosZ.clear(); + fSecondaryNeutronsShieldingIsCaptured.clear(); + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.clear(); + fSecondaryNeutronsShieldingProductionE.clear(); + fSecondaryNeutronsShieldingExitE.clear(); +} +/////////////////////////////////////////////// +/// \brief The main processing event function +/// +TRestEvent* TRestGeant4CosmicNeutronTaggingAnalysisProcess::ProcessEvent(TRestEvent* evInput) { + fInputG4Event = (TRestGeant4Event*)evInput; + *fOutputG4Event = *((TRestGeant4Event*)evInput); + + Reset(); + std::map volume_energy_map; + + for (unsigned int i = 0; i < fVetoVolumeIds.size(); i++) { + int id = fVetoVolumeIds[i]; + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + + Double_t energy = fOutputG4Event->GetEnergyDepositedInVolume(id); + volume_energy_map[volume_name] = energy; + } + + Double_t energy_veto_max = 0; + for (const auto& pair : volume_energy_map) { + Double_t veto_energy = pair.second; + SetObservableValue(pair.first + "VolumeEDep", veto_energy); + if (veto_energy > energy_veto_max) { + energy_veto_max = veto_energy; + }; + } + SetObservableValue("vetoAllEVetoMax", energy_veto_max); + + // veto groups + for (const auto& pair : fVetoGroupVolumeNames) { + Double_t energy_veto_max_group = 0; + for (unsigned int i = 0; i < pair.second.size(); i++) { + string volume_name = pair.second[i]; + Double_t veto_energy = volume_energy_map[volume_name]; + if (veto_energy > energy_veto_max_group) { + energy_veto_max_group = veto_energy; + }; + } + // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) + string group_name; + for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { + if (it == pair.first.cbegin()) { + group_name += std::toupper(*it); + } else { + group_name += std::tolower(*it); + } + } + SetObservableValue("vetoGroup" + group_name + "EVetoMax", energy_veto_max_group); + } + + // + + for (const auto& quenching_factor : fQuenchingFactors) { + string quenching_factor_string = std::to_string(quenching_factor); + // replace "." in string by "_" because its gives very strange problems + quenching_factor_string = + quenching_factor_string.replace(quenching_factor_string.find("."), sizeof(".") - 1, "_"); + volume_energy_map.clear(); + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + for (const auto& id : fVetoVolumeIds) { + string volume_name = (string)fG4Metadata->GetActiveVolumeName(id); + + if (particle_name == "e-" || particle_name == "e+" || particle_name == "gamma") { + // no quenching factor + volume_energy_map[volume_name] += track->GetEnergyInVolume(id); + } else { + // apply quenching factor + volume_energy_map[volume_name] += quenching_factor * track->GetEnergyInVolume(id); + } + } + } + + Double_t energy_veto_max = 0; + for (const auto& pair : volume_energy_map) { + Double_t veto_energy = pair.second; + SetObservableValue(pair.first + "VolumeEDep" + "Qf" + quenching_factor_string, veto_energy); + if (veto_energy > energy_veto_max) { + energy_veto_max = veto_energy; + }; + } + SetObservableValue(string("vetoAllEVetoMax") + "Qf" + quenching_factor_string, energy_veto_max); + + // veto groups + for (const auto& pair : fVetoGroupVolumeNames) { + Double_t energy_veto_max_group = 0; + for (unsigned int i = 0; i < pair.second.size(); i++) { + string volume_name = pair.second[i]; + Double_t veto_energy = volume_energy_map[volume_name]; + if (veto_energy > energy_veto_max_group) { + energy_veto_max_group = veto_energy; + }; + } + // convert to Upper + lower case (VetoGroupTopEVetoMax, ...) + string group_name; + for (auto it = pair.first.cbegin(); it != pair.first.cend(); ++it) { + if (it == pair.first.cbegin()) { + group_name += std::toupper(*it); + } else { + group_name += std::tolower(*it); + } + } + SetObservableValue("vetoGroup" + group_name + "EVetoMax" + "Qf" + quenching_factor_string, + energy_veto_max_group); + } + } + + std::set neutronsCaptured = {}; + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "neutron") { + auto hits = track->GetHits(); + for (int j = 0; j < hits->GetNumberOfHits(); j++) { + string process_name = (string)track->GetProcessName(hits->GetProcess(j)); + if (process_name == "nCapture") { + // << "Neutron capture!!!!!! " << particle_name << "trackId " << track->GetTrackID() + // << " hit " << j << endl; + // track->PrintTrack(); + // hits->PrintHits(j + 1); + + neutronsCaptured.insert(track->GetTrackID()); + + fNeutronsCapturedNumber += 1; + fNeutronsCapturedPosX.push_back(hits->GetX(j)); + fNeutronsCapturedPosY.push_back(hits->GetY(j)); + fNeutronsCapturedPosZ.push_back(hits->GetZ(j)); + + Int_t volumeId = hits->GetVolumeId(j); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + fNeutronsCapturedIsCaptureVolume.push_back(isCaptureVolume); + fNeutronsCapturedProductionE.push_back(track->GetKineticEnergy()); + + // get energy deposited by neutron that undergoes capture and children + double neutronsCapturedEDepByNeutron = 0; + double neutronsCapturedEDepByNeutronAndChildren = 0; + double neutronsCapturedEDepByNeutronInVeto = 0; + double neutronsCapturedEDepByNeutronAndChildrenInVeto = 0; + + std::set parents = {track->GetTrackID()}; + std::map energy_in_veto; + for (int child = 0; child < fOutputG4Event->GetNumberOfTracks(); child++) { + auto track_child = fOutputG4Event->GetTrack(child); + if ((parents.count(track_child->GetParentID()) > 0) || + parents.count(track_child->GetTrackID()) > 0) { + // track or parent is in list of tracks, we add to list and add energy + parents.insert(track_child->GetTrackID()); + neutronsCapturedEDepByNeutronAndChildren += track_child->GetEnergy(); + if (track_child->GetTrackID() == track->GetTrackID()) { + neutronsCapturedEDepByNeutron += track_child->GetEnergy(); + } + for (const auto& vetoId : fVetoVolumeIds) { + neutronsCapturedEDepByNeutronAndChildrenInVeto += + track_child->GetEnergyInVolume(vetoId); + energy_in_veto[vetoId] += track_child->GetEnergyInVolume(vetoId); + if (track_child->GetTrackID() == track->GetTrackID()) { + neutronsCapturedEDepByNeutronInVeto += + track_child->GetEnergyInVolume(vetoId); + } + } + } + } + + fNeutronsCapturedEDepByNeutron.push_back(neutronsCapturedEDepByNeutron); + fNeutronsCapturedEDepByNeutronAndChildren.push_back( + neutronsCapturedEDepByNeutronAndChildren); + fNeutronsCapturedEDepByNeutronInVeto.push_back(neutronsCapturedEDepByNeutronInVeto); + fNeutronsCapturedEDepByNeutronAndChildrenInVeto.push_back( + neutronsCapturedEDepByNeutronAndChildrenInVeto); + + // get max and min energy in each veto (to compare with energy in ALL vetoes) + double energyMaxVeto = 0; + double energyMinVeto = -1; + for (const auto& pair : energy_in_veto) { + auto E = pair.second; + if (E > energyMaxVeto) energyMaxVeto = E; + if (E < energyMaxVeto || energyMinVeto == -1) energyMinVeto = E; + } + + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax.push_back(energyMaxVeto); + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin.push_back(energyMinVeto); + } + } + } + } + + SetObservableValue("neutronsCapturedNumber", fNeutronsCapturedNumber); + SetObservableValue("neutronsCapturedPosX", fNeutronsCapturedPosX); + SetObservableValue("neutronsCapturedPosY", fNeutronsCapturedPosY); + SetObservableValue("neutronsCapturedPosZ", fNeutronsCapturedPosZ); + SetObservableValue("neutronsCapturedIsCaptureVolume", fNeutronsCapturedIsCaptureVolume); + SetObservableValue("neutronsCapturedProductionE", fNeutronsCapturedProductionE); + SetObservableValue("neutronsCapturedEDepByNeutron", fNeutronsCapturedEDepByNeutron); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildren", fNeutronsCapturedEDepByNeutronAndChildren); + SetObservableValue("neutronsCapturedEDepByNeutronInVeto", fNeutronsCapturedEDepByNeutronInVeto); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVeto", + fNeutronsCapturedEDepByNeutronAndChildrenInVeto); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMax", + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMax); + SetObservableValue("neutronsCapturedEDepByNeutronAndChildrenInVetoMin", + fNeutronsCapturedEDepByNeutronAndChildrenInVetoMin); + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "gamma") { + // check if gamma is child of captured neutron + Int_t parent = track->GetParentID(); + if (neutronsCaptured.count(parent) > 0) { + auto hits = track->GetHits(); + + fGammasNeutronCaptureNumber += 1; + fGammasNeutronCapturePosX.push_back(hits->GetX(0)); + fGammasNeutronCapturePosY.push_back(hits->GetY(0)); + fGammasNeutronCapturePosZ.push_back(hits->GetZ(0)); + + Int_t volumeId = hits->GetVolumeId(0); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + fGammasNeutronCaptureIsCaptureVolume.push_back(isCaptureVolume); + fGammasNeutronCaptureProductionE.push_back(track->GetKineticEnergy()); + + // cout << "gamma capture" << endl; + + // hits->PrintHits(1); + } + } + } + + SetObservableValue("gammasNeutronCaptureNumber", fGammasNeutronCaptureNumber); + SetObservableValue("gammasNeutronCapturePosX", fGammasNeutronCapturePosX); + SetObservableValue("gammasNeutronCapturePosY", fGammasNeutronCapturePosY); + SetObservableValue("gammasNeutronCapturePosZ", fGammasNeutronCapturePosZ); + SetObservableValue("gammasNeutronCaptureIsCaptureVolume", fGammasNeutronCaptureIsCaptureVolume); + SetObservableValue("gammasNeutronCaptureProductionE", fGammasNeutronCaptureProductionE); + + std::set secondaryNeutrons = {}; // avoid counting twice + for (int i = 0; i < fOutputG4Event->GetNumberOfTracks(); i++) { + auto track = fOutputG4Event->GetTrack(i); + string particle_name = (string)track->GetParticleName(); + if (particle_name == "neutron" && track->GetParentID() != 0) { // not consider primary + // check if neutron exits shielding + auto hits = track->GetHits(); + for (int j = 0; j < hits->GetNumberOfHits(); j++) { + string process_name = (string)track->GetProcessName(hits->GetProcess(j)); + if (process_name == "Transportation") { + for (const auto& id : fShieldingVolumeIds) { + if (hits->GetVolumeId(j) == id) { + // transportation and shielding == exits shielding + if (secondaryNeutrons.count(track->GetTrackID()) == 0) { + // first time adding this secondary neutron + secondaryNeutrons.insert(track->GetTrackID()); + } else { + continue; + } + fSecondaryNeutronsShieldingNumber += 1; + fSecondaryNeutronsShieldingExitPosX.push_back(hits->GetX(j)); + fSecondaryNeutronsShieldingExitPosY.push_back(hits->GetY(j)); + fSecondaryNeutronsShieldingExitPosZ.push_back(hits->GetZ(j)); + + Int_t volumeId = hits->GetVolumeId(j); + Int_t isCaptureVolume = 0; + for (const auto& id : fCaptureVolumeIds) { + if (volumeId == id) { + isCaptureVolume = 1; + continue; + } + } + Int_t isCaptured = 0; + if (neutronsCaptured.count(track->GetTrackID()) > 0) { + isCaptured = 1; + } + fSecondaryNeutronsShieldingIsCaptured.push_back(isCaptured); + if (isCaptured) + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back( + isCaptureVolume); + else { + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume.push_back(0); + } + + fSecondaryNeutronsShieldingProductionE.push_back(track->GetKineticEnergy()); + fSecondaryNeutronsShieldingExitE.push_back(hits->GetKineticEnergy(j)); + } + } + } + } + } + } + + SetObservableValue("secondaryNeutronsShieldingNumber", fSecondaryNeutronsShieldingNumber); + SetObservableValue("secondaryNeutronsShieldingExitPosX", fSecondaryNeutronsShieldingExitPosX); + SetObservableValue("secondaryNeutronsShieldingExitPosY", fSecondaryNeutronsShieldingExitPosY); + SetObservableValue("secondaryNeutronsShieldingExitPosZ", fSecondaryNeutronsShieldingExitPosZ); + SetObservableValue("secondaryNeutronsShieldingIsCaptured", fSecondaryNeutronsShieldingIsCaptured); + SetObservableValue("secondaryNeutronsShieldingIsCapturedInCaptureVolume", + fSecondaryNeutronsShieldingIsCapturedInCaptureVolume); + SetObservableValue("secondaryNeutronsShieldingProductionE", fSecondaryNeutronsShieldingProductionE); + SetObservableValue("secondaryNeutronsShieldingExitE", fSecondaryNeutronsShieldingExitE); + + return fOutputG4Event; +} + +/////////////////////////////////////////////// +/// \brief Function to include required actions after all events have been +/// processed. +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::EndProcess() { + // Function to be executed once at the end of the process + // (after all events have been processed) + + // Start by calling the EndProcess function of the abstract class. + // Comment this if you don't want it. + // TRestEventProcess::EndProcess(); +} + +/////////////////////////////////////////////// +/// \brief Function to read input parameters from the RML +/// TRestGeant4CosmicNeutronTaggingAnalysisProcess metadata section +/// +void TRestGeant4CosmicNeutronTaggingAnalysisProcess::InitFromConfigFile() { + // word to identify active volume as veto (default = "veto" e.g. "vetoTop") + string veto_keyword = GetParameter("vetoKeyword", "veto"); + fVetoKeyword = clean_string(veto_keyword); + // comma separated tags: "top, bottom, ..." + string veto_group_keywords = GetParameter("vetoGroupKeywords", ""); + stringstream ss(veto_group_keywords); + while (ss.good()) { + string substr; + getline(ss, substr, ','); + fVetoGroupKeywords.push_back(clean_string(substr)); + } + + // word to identify active volume as capture sheet (cadmium, default = "sheet" e.g. + // "scintillatorSheetTop1of4") + string capture_keyword = GetParameter("captureKeyword", "sheet"); + fCaptureKeyword = clean_string(capture_keyword); + + // word to identify active volume as shielding + + string shielding_keyword = GetParameter("shieldingKeyword", "shielding"); + fShieldingKeyword = clean_string(shielding_keyword); + + // comma separated quenching factors: "0.15, 1.00, ..." + string quenching_factors = GetParameter("vetoQuenchingFactors", "-1"); + stringstream ss_qf(quenching_factors); + while (ss_qf.good()) { + string substr; + getline(ss_qf, substr, ','); + substr = clean_string(substr); + Float_t quenching_factor = (Float_t)std::atof(substr.c_str()); + if (quenching_factor > 1 || quenching_factor < 0) { + cout << "ERROR: quenching factor must be between 0 and 1" << endl; + continue; + } + fQuenchingFactors.push_back(quenching_factor); + } +} diff --git a/src/TRestGeant4VetoAnalysisProcess.cxx b/src/TRestGeant4VetoAnalysisProcess.cxx index 2a28410c..3c601192 100644 --- a/src/TRestGeant4VetoAnalysisProcess.cxx +++ b/src/TRestGeant4VetoAnalysisProcess.cxx @@ -21,8 +21,34 @@ *************************************************************************/ ////////////////////////////////////////////////////////////////////////// -/// TODO: fill documentation +/// The TRestGeant4VetoAnalysisProcess generates observables based on veto volumes energy depositions. It was +/// first developed as a process for the IAXO experiment but can be used in any analysis. It uses `keywords` +/// to identify different relevant volumes (such as vetoes). The default veto keyword for IAXO is `veto` and +/// it will tag each volume containing the keyword as a veto volume, so avoid using the keyword on volumes +/// that do not act as vetoes. /// +/// Parameters: +/// `vetoKeyword`: keyword to identify a volume of the geometry as a veto. The condition is that the keyword +/// is contained inside the name of the volume. Only volumes serving as veto (i.e. scintillators) should +/// contain this vetoKeyword in their names. +/// `vetoGroupKeywords`: comma separated keywords used to identify different groups of vetoes. This is an +/// optional parameter that when specified will make the process return additional observables on a per group +/// basis. The most common use case is using group names to identify the location of the vetoes (top, bottom, +/// ...). The volumes detected as vetoes (via vetoKeyword) will also be assigned to a group if they contain a +/// keyword contained in this list. +/// `vetoQuenchingFactors`: comma separated values for the quenching factors used in the analysis. The +/// observables will be calculated for each of the quenching factors contained in this list. Values between 0 +/// and 1 only. This is useful in case the user doesn't know the exact value of the quenching factor. Also it +/// is useful to insert `0` or `1` to study the effects of electromagnetic processes only, or no quenching at +/// all. +/// +/// Example usage: +/// +/// +/// +/// +/// ///-------------------------------------------------------------------------- /// /// RESTsoft - Software for Rare Event Searches with TPCs