diff --git a/CMakeCache.txt b/CMakeCache.txt deleted file mode 100644 index c07fc400..00000000 --- a/CMakeCache.txt +++ /dev/null @@ -1,38 +0,0 @@ -# This is the CMakeCache file. -# For build in directory: /tmp/locust_mc -# It was generated by CMake: /usr/bin/cmake -# You can edit this file to change values found and used by cmake. -# If you do not want to change any of the values, simply exit the editor. -# If you do want to change a value, simply edit, save, and exit the editor. -# The syntax for the file is as follows: -# KEY:TYPE=VALUE -# KEY is the name of a variable in the cache. -# TYPE is a hint to GUIs for the type of VALUE, DO NOT EDIT TYPE!. -# VALUE is the current value for the KEY. - -######################## -# EXTERNAL cache entries -######################## - - -######################## -# INTERNAL cache entries -######################## - -//This is the directory where this CMakeCache.txt was created -CMAKE_CACHEFILE_DIR:INTERNAL=/tmp/locust_mc -//Major version of cmake used to create the current loaded cache -CMAKE_CACHE_MAJOR_VERSION:INTERNAL=3 -//Minor version of cmake used to create the current loaded cache -CMAKE_CACHE_MINOR_VERSION:INTERNAL=22 -//Patch version of cmake used to create the current loaded cache -CMAKE_CACHE_PATCH_VERSION:INTERNAL=1 -//Path to CMake executable. -CMAKE_COMMAND:INTERNAL=/usr/bin/cmake -//Path to cpack program executable. -CMAKE_CPACK_COMMAND:INTERNAL=/usr/bin/cpack -//Path to ctest program executable. -CMAKE_CTEST_COMMAND:INTERNAL=/usr/bin/ctest -//Path to CMake installation. -CMAKE_ROOT:INTERNAL=/usr/share/cmake-3.22 - diff --git a/CMakeLists.txt b/CMakeLists.txt index 8470f4f8..684a9b3e 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -6,7 +6,7 @@ cmake_minimum_required( VERSION 3.1 ) # Define the project cmake_policy( SET CMP0048 NEW ) # version in project() -project( locust_mc VERSION 3.0.0) +project( locust_mc VERSION 3.1.0) list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/Scarab/cmake ) @@ -189,12 +189,18 @@ if (ROOT_FOUND) add_definitions(-DROOT_FOUND) pbuilder_add_ext_libraries (${ROOT_LIBRARIES}) message(STATUS "ROOT_LIBRARIES is ${ROOT_LIBRARIES}") + message(STATUS "ROOT_INCLUDE_DIRS is ${ROOT_INCLUDE_DIRS}") else (ROOT_FOUND) message(STATUS "Building without ROOT") remove_definitions(-DROOT_FOUND) endif (ROOT_FOUND) include_directories (${ROOT_INCLUDE_DIR}) +option (locust_mc_USE_MPI "Use MPI to accelerate calculations using multiple processors" OFF) +if (locust_mc_USE_MPI) + find_package(MPI REQUIRED) +endif (locust_mc_USE_MPI) + # Boost (1.46 required for filesystem version 3) list (APPEND Boost_COMPONENTS date_time filesystem program_options system thread) find_package (Boost 1.46.0 REQUIRED COMPONENTS ${Boost_COMPONENTS}) diff --git a/Source/Applications/Testing/testFFT.cc b/Source/Applications/Testing/testFFT.cc index 5d45606d..7140255d 100644 --- a/Source/Applications/Testing/testFFT.cc +++ b/Source/Applications/Testing/testFFT.cc @@ -73,6 +73,10 @@ double GetPower(bool bTime) ptransform += data[j][0]*data[j][0]+data[j][1]*data[j][1]; } + fftw_destroy_plan( plan); + fftw_free( data ); + data = NULL; + if (bTime) { LPROG(testlog, "power of original data time series is: " << pdata); diff --git a/Source/Applications/Testing/testLMCTestSignal.cc b/Source/Applications/Testing/testLMCTestSignal.cc index 82966b80..14354d08 100644 --- a/Source/Applications/Testing/testLMCTestSignal.cc +++ b/Source/Applications/Testing/testLMCTestSignal.cc @@ -113,6 +113,10 @@ class testLMCTestSignal ptransform += data[j][0]*data[j][0]+data[j][1]*data[j][1]; } + fftw_destroy_plan( plan); + fftw_free( data ); + data = NULL; + if (bTime) { LPROG(testlog, "fft power of original data time series is: " << pdata); @@ -161,6 +165,7 @@ TEST_CASE( "LMCTestSignal with default parameter values (pass)", "[single-file]" double threshold = 1.e-4; REQUIRE( fabs( aTestLMCTestSignal.GetPower(aSignal, N0, 1) - aTestLMCTestSignal.GetPower(aSignal, N0, 0) ) < threshold*aTestLMCTestSignal.GetPower(aSignal, N0, 1) ); REQUIRE( fabs( aTestLMCTestSignal.GetPower(aSignal, N0, 1)/N0/50. - pow(aTestLMCTestSignal.GetAmplitude(),2.) ) < threshold*aTestLMCTestSignal.GetPower(aSignal, N0, 1)/N0/50.); + aSignal->Reset(); } diff --git a/Source/Applications/Testing/testMockEGun.cc b/Source/Applications/Testing/testMockEGun.cc index 3d657205..c8b153fb 100644 --- a/Source/Applications/Testing/testMockEGun.cc +++ b/Source/Applications/Testing/testMockEGun.cc @@ -103,7 +103,10 @@ double GetPower() } LPROG(testlog, "E-gun transformed data sum is: " << ptransform/N0); + fftw_destroy_plan(plan); + fftw_free( data ); + data = NULL; return ptransform/N0; diff --git a/Source/Applications/Testing/testMockFreeField.cc b/Source/Applications/Testing/testMockFreeField.cc index d7590bd3..5b38fb7a 100644 --- a/Source/Applications/Testing/testMockFreeField.cc +++ b/Source/Applications/Testing/testMockFreeField.cc @@ -164,6 +164,8 @@ TEST_CASE( "Mock free space Larmor power. (pass)", "[single-file]" ) fftw_destroy_plan(plan); + fftw_free( data ); + data = NULL; double threshold = 1.e-4; diff --git a/Source/CMakeLists.txt b/Source/CMakeLists.txt index 6209d4b6..95753f98 100644 --- a/Source/CMakeLists.txt +++ b/Source/CMakeLists.txt @@ -13,6 +13,7 @@ set( LOCUST_MC_HEADER_FILES Core/LMCHFSSResponseFileHandler.hh Core/LMCFIRFileHandler.hh Core/LMCTFFileHandler.hh + Core/LMCMPIEnvironment.hh Utilities/LMCUtility.hh Utilities/LMCCavityUtility.hh @@ -34,8 +35,6 @@ set( LOCUST_MC_HEADER_FILES Generators/LMCTestSignalGenerator.hh Generators/LMCDecimateSignalGenerator.hh Generators/LMCLowPassFilterFFTGenerator.hh - Generators/LMCHighPassFilterFFTGenerator.hh - Generators/LMCLocalOscillatorGenerator.hh Generators/LMCButterworthLPFGenerator.hh Transmitters/LMCFieldBuffer.hh @@ -108,8 +107,6 @@ set( LOCUST_MC_SOURCE_FILES Generators/LMCTestSignalGenerator.cc Generators/LMCDecimateSignalGenerator.cc Generators/LMCLowPassFilterFFTGenerator.cc - Generators/LMCHighPassFilterFFTGenerator.cc - Generators/LMCLocalOscillatorGenerator.cc Generators/LMCButterworthLPFGenerator.cc @@ -146,16 +143,25 @@ set( LOCUST_MC_SOURCE_FILES ) +if (locust_mc_USE_MPI) + set( LOCUST_MC_HEADER_FILES ${LOCUST_MC_HEADER_FILES} + Core/LMCMPIInterface.hh + ) + + set( LOCUST_MC_SOURCE_FILES ${LOCUST_MC_SOURCE_FILES} + Core/LMCMPIInterface.cc + ) +endif (locust_mc_USE_MPI) + + if (locust_mc_BUILD_WITH_KASSIOPEIA) # Here we only need to add the Kassiopeia-dependent header/source files to the header/source-file lists set( LOCUST_MC_HEADER_FILES ${LOCUST_MC_HEADER_FILES} - Generators/LMCFreeFieldSignalGenerator.hh Generators/LMCKassSignalGenerator.hh Generators/LMCArraySignalGenerator.hh - Generators/LMCCavitySignalGenerator.hh - + Generators/LMCCavitySignalGenerator.hh Transmitters/LMCKassTransmitter.hh Transmitters/LMCKassCurrentTransmitter.hh diff --git a/Source/Core/LMCFIRFileHandler.cc b/Source/Core/LMCFIRFileHandler.cc index e7e6fd4a..d8b1f322 100644 --- a/Source/Core/LMCFIRFileHandler.cc +++ b/Source/Core/LMCFIRFileHandler.cc @@ -24,10 +24,10 @@ namespace locust bool FIRTransmitterHandler::Configure(const scarab::param_node& aParam) { - int fNModes = 2; - if( aParam.has( "n-modes" ) ) + int fNModes = 2; + if( aParam.has( "n-modes" ) ) { - fNModes = aParam["n-modes"]().as_int(); + fNModes = aParam["n-modes"]().as_int(); } if( aParam.has( "fir-transmitter-filename" ) ) { diff --git a/Source/Core/LMCHFSSResponseFileHandler.cc b/Source/Core/LMCHFSSResponseFileHandler.cc index 03169e31..62f002aa 100644 --- a/Source/Core/LMCHFSSResponseFileHandler.cc +++ b/Source/Core/LMCHFSSResponseFileHandler.cc @@ -175,7 +175,24 @@ namespace locust fftw_free(fFIRComplex); fFIRComplex = NULL; } + + + for( int bTE=0; bTE<2; bTE++) + { + for(int l=0; lGetProcess() == 0) +#define MPI_SECOND_PROCESS if (KEMField::KMPIInterface::GetInstance()->GetProcess() == 1) +#else +#define MPI_SINGLE_PROCESS +#define MPI_SECOND_PROCESS +#endif + +#endif /* KMPIENVIRONMENT_HH_ */ diff --git a/Source/Core/LMCMPIInterface.cc b/Source/Core/LMCMPIInterface.cc new file mode 100644 index 00000000..d0be74b0 --- /dev/null +++ b/Source/Core/LMCMPIInterface.cc @@ -0,0 +1,428 @@ +#include "LMCMPIInterface.hh" + +//#include "KEMCoreMessage.hh" + +#ifdef LOCAL_RANK_MPI +//used to retrieve the host name +//should be available on on POSIX systems +#include +#include +#endif + +#define MSG_TAG 999 +#define HOST_DETERMINATION_TAG 998 +#define LOCALID_DETERMINATION_TAG 997 + +namespace locust +{ +LMCMPIInterface* LMCMPIInterface::fMPIInterface = nullptr; + +LMCMPIInterface::LMCMPIInterface() +{ + fProcess = -1; + fNProcesses = -1; + fLocalRank = -1; + fSubGroupRank = -1; + fNSubGroupProcesses = -1; + fPartnerProcessID = -1; + fSplitMode = false; +} + +LMCMPIInterface::~LMCMPIInterface() = default; + +void LMCMPIInterface::SetSplitMode(bool enable) +{ + fSplitMode = enable; + + //construct groups/communicators to evenly split up processes + SetupSubGroups(); + + //cannot make a valid split, so revert to standard mode + if (fSplitMode && !fValidSplit) { +// kem_cout(eWarning) << "Invalid MPI split mode detected - revert to standard mode with " << fNProcesses << " processes." << eom; + fSplitMode = false; + } + else if (fSplitMode) { +// kem_cout(eNormal) << "Using MPI split mode with 2 x " << fNSubGroupProcesses << " processes" << eom; + } +} + +void LMCMPIInterface::Initialize(int* argc, char*** argv, bool split_mode) +{ + /* Let the system do what it needs to start up MPI */ + int initialized = 0; + MPI_Initialized(&initialized); + + if (!initialized) + MPI_Init(argc, argv); + + /* Get my process fProcess */ + MPI_Comm_rank(MPI_COMM_WORLD, &fProcess); + + /* Find out how many processes are being used */ + MPI_Comm_size(MPI_COMM_WORLD, &fNProcesses); + + if (fProcess <= 0) { + if (fNProcesses <= 0) { +// kem_cout(eWarning) << "No MPI processes found - not running in an MPI context?" << eom; + } + else if (!initialized) { // only show this once +// kem_cout << "Running in MPI context with " << fNProcesses << " processes." << eom; + } + } + + //now determine the local rank of this process (indexed from zero) on the local host + //for example, if processes (0,2,5) are running on host A + //and processes (1,3,4) are running on host B, then the + //local rank of process 3 is 1, and the local rank of process 5 is 2 + DetermineLocalRank(); + + SetSplitMode(split_mode); +} + +void LMCMPIInterface::Finalize() +{ + /* Shut down MPI */ + int finalized = 0; + MPI_Finalized(&finalized); + + if (!finalized) + MPI_Finalize(); +} + +/** + * Interface to accessing LMCMPIInterface. + */ +LMCMPIInterface* LMCMPIInterface::GetInstance() +{ + if (fMPIInterface == nullptr) + fMPIInterface = new LMCMPIInterface(); + return fMPIInterface; +} + +/** + * Ensures that a process written between BeginSequentialProcess() and + * EndSequentialProcess() is done one processor at a time. + */ +void LMCMPIInterface::BeginSequentialProcess() +{ + int flag = 1; + + GlobalBarrier(); + + if (fProcess > 0) + MPI_Recv(&flag, 1, MPI_INT, fProcess - 1, 50, MPI_COMM_WORLD, &fStatus); +} + +/** + * @see BeginSequentialProcess() + */ +void LMCMPIInterface::EndSequentialProcess() +{ + int flag = 0; + + if (fProcess < (fNProcesses - 1)) + MPI_Send(&flag, 1, MPI_INT, fProcess + 1, 50, MPI_COMM_WORLD); + + GlobalBarrier(); +} + +//safely print messages in process order, by passing information to the +//root process for collection before calling cout +void LMCMPIInterface::PrintMessage(std::string msg) +{ + if (fNProcesses < 0 || fProcess < 0) { +// kem_cout(eError) << "MPI not initialized" << ret; +// kem_cout << "Message was: " << msg << eom; + return; + } + + unsigned int n_char = msg.size(); + + std::vector in_msg_sizes; + std::vector out_msg_sizes; + in_msg_sizes.resize(fNProcesses, 0); + out_msg_sizes.resize(fNProcesses, 0); + in_msg_sizes[fProcess] = n_char; + + //obtain the message sizes from all of the objects + MPI_Allreduce(&(in_msg_sizes[0]), &(out_msg_sizes[0]), fNProcesses, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + + //compute the total message size + unsigned int total_msg_size = 0; + std::vector msg_start_indexes; + msg_start_indexes.resize(fNProcesses); + for (int i = 0; i < fNProcesses; i++) { + total_msg_size += out_msg_sizes[i]; + }; + + for (int i = 0; i < fNProcesses; i++) { + for (int j = 0; j < i; j++) { + msg_start_indexes[i] += out_msg_sizes[j]; + } + }; + + //allocate buffers to reduce all of the messages + std::vector buf; + buf.resize(total_msg_size); + + //fill the appropriate section of the buffer + for (unsigned int i = 0; i < msg.size(); i++) { + buf[msg_start_indexes[fProcess] + i] = msg.at(i); + } + + MPI_Status status; + if (fProcess == 0) { + for (int i = 1; i < fNProcesses; i++) { + MPI_Recv(&(buf[msg_start_indexes[i]]), out_msg_sizes[i], MPI_CHAR, i, MSG_TAG, MPI_COMM_WORLD, &status); + } + } + else { + MPI_Send(&(buf[msg_start_indexes[fProcess]]), out_msg_sizes[fProcess], MPI_CHAR, 0, MSG_TAG, MPI_COMM_WORLD); + } + + if (fProcess == 0) { + //convert to string + std::stringstream final_output; + for (char c : buf) { + final_output << c; + } + std::string full_message = final_output.str(); + + //split the messages into separate lines so we can correctly use KEMField::cout + std::vector lines; + std::string::size_type pos = 0; + std::string::size_type prev = 0; + while ((pos = full_message.find('\n', prev)) != std::string::npos) { + lines.push_back(full_message.substr(prev, pos - prev)); + prev = pos + 1; + } + + if (lines.size() == 0) { + //message has no newline characters, so push the whole message into one line + lines.push_back(full_message); + } + + //print message line by line + for (auto& line : lines) { +// kem_cout() << line << eom; + } + } +} + + +void LMCMPIInterface::DetermineLocalRank() +{ +#ifdef LOCAL_RANK_MPI + //get the machine's hostname + char host_name[256]; + int ret_val = gethostname(host_name, 256); + if (ret_val != 0) { + std::cout << "Host name error!" << std::endl; + LMCMPIInterface::GetInstance()->Finalize(); + std::exit(1); + }; + + std::stringstream hostname_ss; + int count = 0; + do { + hostname_ss << host_name[count]; + count++; + } while (host_name[count] != '\0' && count < 256); + std::string hostname = hostname_ss.str(); + fHostName = hostname; + + //first we have to collect all of the hostnames that are running a process + unsigned int n_char = hostname.size(); + std::vector in_msg_sizes; + std::vector out_msg_sizes; + in_msg_sizes.resize(fNProcesses, 0); + out_msg_sizes.resize(fNProcesses, 0); + in_msg_sizes[fProcess] = n_char; + + //obtain the message sizes from all of the processes + MPI_Allreduce(&(in_msg_sizes[0]), &(out_msg_sizes[0]), fNProcesses, MPI_UNSIGNED, MPI_SUM, MPI_COMM_WORLD); + + //compute the total message size + unsigned int total_msg_size = 0; + std::vector msg_start_indexes; + msg_start_indexes.resize(fNProcesses); + for (int i = 0; i < fNProcesses; i++) { + total_msg_size += out_msg_sizes[i]; + }; + + for (int i = 0; i < fNProcesses; i++) { + for (int j = 0; j < i; j++) { + msg_start_indexes[i] += out_msg_sizes[j]; + } + }; + + //allocate buffers to reduce all of the messages + std::vector buf; + buf.resize(total_msg_size); + + //fill the appropriate section of the buffer + for (unsigned int i = 0; i < hostname.size(); i++) { + buf[msg_start_indexes[fProcess] + i] = hostname.at(i); + } + + //reduce the buffer across all processes + MPI_Status status; + if (fProcess == 0) { + for (int i = 1; i < fNProcesses; i++) { + MPI_Recv(&(buf[msg_start_indexes[i]]), + out_msg_sizes[i], + MPI_CHAR, + i, + HOST_DETERMINATION_TAG, + MPI_COMM_WORLD, + &status); + } + } + else { + unsigned int root_rank = 0; + MPI_Send(&(buf[msg_start_indexes[fProcess]]), + out_msg_sizes[fProcess], + MPI_CHAR, + root_rank, + HOST_DETERMINATION_TAG, + MPI_COMM_WORLD); + } + + //now broadcast the complete list of hostnames to all processes + MPI_Bcast(&(buf[0]), total_msg_size, MPI_CHAR, 0, MPI_COMM_WORLD); + + //now every node has a list of all host names + //we now need to figure out how many other processes are also running on + //the same host, and how many devices are available on this host + //then we can distribute the processes equitably to each device + + std::vector hostname_list; + hostname_list.resize(fNProcesses); + for (int i = 0; i < fNProcesses; i++) { + hostname_list[i] = std::string(""); + for (unsigned int j = 0; j < out_msg_sizes[i]; j++) { + hostname_list[i].push_back(buf[msg_start_indexes[i] + j]); + } + } + + //collect all the process ids of all the process running on this host + fCoHostedProcessIDs.clear(); + for (int i = 0; i < fNProcesses; i++) { + if (hostname == hostname_list[i]) { + fCoHostedProcessIDs.push_back(i); + } + } + + //determine the 'local' rank of this process + for (unsigned int i = 0; i < fCoHostedProcessIDs.size(); i++) { + if (fCoHostedProcessIDs[i] == fProcess) { + fLocalRank = i; + }; + } + +#endif +} + + +void LMCMPIInterface::SetupSubGroups() +{ +#ifdef LOCAL_RANK_MPI + //we need to retrieve the local rank from each process + //to make a associative map betweek global-rank and local-rank + std::vector local_ranks; + local_ranks.resize(fNProcesses); + local_ranks[fProcess] = fLocalRank; + + //reduce the buffer across all processes + MPI_Status status; + if (fProcess == 0) { + for (int i = 1; i < fNProcesses; i++) { + MPI_Recv(&(local_ranks[i]), 1, MPI_INT, i, LOCALID_DETERMINATION_TAG, MPI_COMM_WORLD, &status); + } + } + else { + unsigned int root_rank = 0; + MPI_Send(&(local_ranks[fProcess]), 1, MPI_INT, root_rank, LOCALID_DETERMINATION_TAG, MPI_COMM_WORLD); + } + + //now broadcast the complete list of hostnames to all processes + MPI_Bcast(&(local_ranks[0]), fNProcesses, MPI_INT, 0, MPI_COMM_WORLD); + + //now every process has a list of the local rank associated with every other process + //now we can proceed to determine which group they below to + std::vector even_members; + std::vector odd_members; + for (int i = 0; i < fNProcesses; i++) { + if (local_ranks[i] % 2 == 0) { + even_members.push_back(i); + } + else { + odd_members.push_back(i); + } + } + + fValidSplit = false; + if (even_members.size() == odd_members.size()) { + fValidSplit = true; + } + + //get the world group + MPI_Group world; + MPI_Comm_group(MPI_COMM_WORLD, &world); + + //now we go ahead and construct the groups and communicators + MPI_Group_incl(world, even_members.size(), &(even_members[0]), &fEvenGroup); + MPI_Group_incl(world, odd_members.size(), &(odd_members[0]), &fOddGroup); + + MPI_Comm_create(MPI_COMM_WORLD, fEvenGroup, &fEvenCommunicator); + MPI_Comm_create(MPI_COMM_WORLD, fOddGroup, &fOddCommunicator); + + //now we set things up for this process + if (fLocalRank % 2 == 0) { + fIsEvenGroupMember = true; + MPI_Comm_rank(fEvenCommunicator, &fSubGroupRank); + fNSubGroupProcesses = even_members.size(); + } + else { + fIsEvenGroupMember = false; + MPI_Comm_rank(fOddCommunicator, &fSubGroupRank); + fNSubGroupProcesses = odd_members.size(); + } + + //finally, if we have a valid split (equal numbers of even and odd process) + //we pair up processes so they can exchange data + if (fValidSplit) //must have even number of processes! + { + int status = 0; + int result = 0; + if (fCoHostedProcessIDs.size() % 2 == 0) { + status = 1; + }; + MPI_Allreduce(&status, &result, 1, MPI_INT, MPI_SUM, MPI_COMM_WORLD); + + //to make things faster we first try to pair up processes + //which share the same node/host + if (result == fNProcesses) { + //we can because each host has an even number of processes + if (fIsEvenGroupMember) { + fPartnerProcessID = fCoHostedProcessIDs[fLocalRank + 1]; + } + else { + fPartnerProcessID = fCoHostedProcessIDs[fLocalRank - 1]; + } + } + else { + //this isn't possible so we have to pair up processes across nodes + if (fIsEvenGroupMember) { + fPartnerProcessID = fProcess + 1; + } + else { + fPartnerProcessID = fProcess - 1; + } + } + } +#endif +} + +} // namespace KEMField diff --git a/Source/Core/LMCMPIInterface.hh b/Source/Core/LMCMPIInterface.hh new file mode 100644 index 00000000..e2e15b52 --- /dev/null +++ b/Source/Core/LMCMPIInterface.hh @@ -0,0 +1,153 @@ +#ifndef LMCMPIINTERFACE_DEF +#define LMCMPIINTERFACE_DEF + +#include "mpi.h" + +#include +#include + +#define LOCAL_RANK_MPI + +namespace locust +{ + +class LMCMPIInterface +{ + public: + static LMCMPIInterface* GetInstance(); + + void Initialize(int* argc, char*** argv, bool split_mode = true); + void Finalize(); + + void SetSplitMode(bool enable = true); + + bool Check() const { + return (fProcess >= 0) && (fNProcesses > 0); + } + int GetProcess() const + { + return fProcess; + } + int GetNProcesses() const + { + return fNProcesses; + } + int GetLocalRank() const + { + return fLocalRank; + } + std::string GetHostName() const + { + return fHostName; + }; + + void BeginSequentialProcess(); + void EndSequentialProcess(); + + void GlobalBarrier() const + { + MPI_Barrier(MPI_COMM_WORLD); + } + + //when called, this function must be encountered by all processes + //or the program will lock up, treat as a global barrier + void PrintMessage(std::string msg); + + //routines to be used by programs which split the processes into two + //groups bases on even/odd local process rank + bool SplitMode() + { + return fSplitMode; + }; + bool IsSplitValid() + { + return fValidSplit; + }; + bool IsEvenGroupMember() + { + return fIsEvenGroupMember; + }; + int GetNSubGroupProcesses() + { + return fNSubGroupProcesses; + } + int GetSubGroupRank() + { + return fSubGroupRank; + }; + int GetPartnerProcessID() + { + return fPartnerProcessID; + }; + + MPI_Group* GetSubGroup() + { + if (fIsEvenGroupMember) { + return &fEvenGroup; + } + else { + return &fOddGroup; + }; + } + + MPI_Comm* GetSubGroupCommunicator() + { + if (fIsEvenGroupMember) { + return &fEvenCommunicator; + } + else { + return &fOddCommunicator; + }; + } + + MPI_Group* GetEvenGroup() + { + return &fEvenGroup; + }; + MPI_Group* GetOddGroup() + { + return &fOddGroup; + }; + MPI_Comm* GetEvenCommunicator() + { + return &fEvenCommunicator; + }; + MPI_Comm* GetOddCommunicator() + { + return &fOddCommunicator; + }; + + + protected: + LMCMPIInterface(); + virtual ~LMCMPIInterface(); + + static LMCMPIInterface* fMPIInterface; + + int fProcess; + int fNProcesses; + int fLocalRank; + std::string fHostName; + std::vector fCoHostedProcessIDs; + + //groups and communicators for splitting processes into + //two sets, based on whether they have even/odd (local) ranks + bool fSplitMode; + MPI_Group fEvenGroup; //even process subgroup + MPI_Group fOddGroup; //odd process subgroup + MPI_Comm fEvenCommunicator; //comm for even group + MPI_Comm fOddCommunicator; //comm for odd group + bool fValidSplit; //true if the size of the subgroups is equal + bool fIsEvenGroupMember; //true if this process is a member of the even subgroup + int fSubGroupRank; //rank of this process in its subgroup + int fNSubGroupProcesses; //number of processes in the subgroup this process belongs to + int fPartnerProcessID; //global rank of partner process in other subgroup + + void DetermineLocalRank(); + void SetupSubGroups(); + + MPI_Status fStatus; +}; +} // namespace locust + +#endif /* LMCMPIINTERFACE_DEF */ diff --git a/Source/Generators/LMCHighPassFilterFFTGenerator.cc b/Source/Deprecated/LMCHighPassFilterFFTGenerator.cc similarity index 100% rename from Source/Generators/LMCHighPassFilterFFTGenerator.cc rename to Source/Deprecated/LMCHighPassFilterFFTGenerator.cc diff --git a/Source/Generators/LMCHighPassFilterFFTGenerator.hh b/Source/Deprecated/LMCHighPassFilterFFTGenerator.hh similarity index 100% rename from Source/Generators/LMCHighPassFilterFFTGenerator.hh rename to Source/Deprecated/LMCHighPassFilterFFTGenerator.hh diff --git a/Source/Generators/LMCLocalOscillatorGenerator.cc b/Source/Deprecated/LMCLocalOscillatorGenerator.cc similarity index 100% rename from Source/Generators/LMCLocalOscillatorGenerator.cc rename to Source/Deprecated/LMCLocalOscillatorGenerator.cc diff --git a/Source/Generators/LMCLocalOscillatorGenerator.hh b/Source/Deprecated/LMCLocalOscillatorGenerator.hh similarity index 100% rename from Source/Generators/LMCLocalOscillatorGenerator.hh rename to Source/Deprecated/LMCLocalOscillatorGenerator.hh diff --git a/Source/Fields/LMCCylindricalCavity.cc b/Source/Fields/LMCCylindricalCavity.cc index 7f1864d4..c03c052a 100644 --- a/Source/Fields/LMCCylindricalCavity.cc +++ b/Source/Fields/LMCCylindricalCavity.cc @@ -17,7 +17,9 @@ namespace locust fProbeGain( {1., 1., 1.}), fCavityProbeZ( {0., 0., 0.} ), fCavityProbeRFrac( {0.5, 0.5, 0.5} ), - fCavityProbeTheta( {0.0, 0.0, 0.0} ) + fCavityProbeTheta( {0.0, 0.0, 0.0} ), + fCaterpillarCavity( false ), + fApplyDopplerShift( true ) {} CylindricalCavity::~CylindricalCavity() {} @@ -33,6 +35,16 @@ namespace locust return false; } + if( aParam.has( "caterpillar-cavity" ) ) + { + fCaterpillarCavity = aParam["caterpillar-cavity"]().as_bool(); + } + + if( aParam.has( "apply-doppler-shift" ) ) + { + fApplyDopplerShift = aParam["apply-doppler-shift"]().as_bool(); + } + if( aParam.has( "cavity-radius" ) ) { SetDimR( aParam["cavity-radius"]().as_double() ); @@ -139,6 +151,9 @@ namespace locust SetNormFactors(SetUnityNormFactors(GetNModes(), 0)); // Temporary quick normalization factors of 1.0 } + fFieldCore->ReadBesselZeroes((dataDir / "BesselZeros.txt").string(), 0 ); + fFieldCore->ReadBesselZeroes((dataDir / "BesselPrimeZeros.txt").string(), 1 ); + if( PlotModeMaps() ) { double zSlice = 0.0; @@ -149,8 +164,6 @@ namespace locust PrintModeMaps(GetNModes(), zSlice, thetaSlice); } - fFieldCore->ReadBesselZeroes((dataDir / "BesselZeros.txt").string(), 0 ); - fFieldCore->ReadBesselZeroes((dataDir / "BesselPrimeZeros.txt").string(), 1 ); SetNormFactors(CalculateNormFactors(GetNModes(), 0)); // Calculate the realistic normalization factors. CheckNormalization(GetNModes(), 0); // E fields integration over volume @@ -227,11 +240,16 @@ namespace locust double vz = tKassParticleXP[5]; double term1 = fFieldCore->GetBesselNKPrimeZeros(l,m) / GetDimR(); double term2 = n * LMCConst::Pi() / GetDimL(); + if ( fCaterpillarCavity ) + { + // Assume n-channels is the same as the number of etalon sections: + term2 *= GetNChannels(); + } double lambda = 1. / pow( 1. / 4. / LMCConst::Pi() / LMCConst::Pi() * ( term1*term1 + term2*term2 ), 0.5); double lambda_c = 2 * LMCConst::Pi() * GetDimR() / fFieldCore->GetBesselNKPrimeZeros(l,m); double vp = LMCConst::C() / pow( 1. - lambda*lambda/lambda_c/lambda_c, 0.5 ); double dopplerShift = 0.; - if (vp > 0.) dopplerShift = vz / vp; + if ((vp > 0.) && (fApplyDopplerShift)) dopplerShift = vz / vp; freqPrime.push_back( ( 1. + dopplerShift ) * tKassParticleXP[7] ); return freqPrime; } @@ -325,7 +343,23 @@ namespace locust std::vector tEFieldAtProbe; for (unsigned index=0; index fCavityProbeRFrac; std::vector fCavityProbeTheta; std::vector fProbeGain; + bool fCaterpillarCavity; + bool fApplyDopplerShift; }; diff --git a/Source/Fields/LMCField.cc b/Source/Fields/LMCField.cc index 8fab166c..420ef717 100644 --- a/Source/Fields/LMCField.cc +++ b/Source/Fields/LMCField.cc @@ -16,8 +16,9 @@ namespace locust fCentralFrequency(1.63e11), fAvgDotProductFactor( 0. ), fNModes( 2 ), + fNChannels( 1 ), fbMultiMode( false ), - fTM111( false ), + fTM111( false ), fR( 0.18 ), fL( 3.0 ), fX( 0.010668 ), diff --git a/Source/Fields/LMCModeMapCavity.cc b/Source/Fields/LMCModeMapCavity.cc index 3ab5455d..8a0460ef 100644 --- a/Source/Fields/LMCModeMapCavity.cc +++ b/Source/Fields/LMCModeMapCavity.cc @@ -20,6 +20,7 @@ namespace locust fDim2_max(6.284), fDim3_min(0.), fDim3_max(0.1524), + fZshift(0.), fnPixel1(10), fnPixel2(10), fnPixel3(10) @@ -73,6 +74,11 @@ namespace locust { fDim3_max = aParam["dim3-max"]().as_double(); } + if( aParam.has( "mode-map-z-shift" ) ) + { + fZshift = aParam["mode-map-z-shift"]().as_double(); + } + return true; } @@ -116,7 +122,7 @@ namespace locust std::string token; std::stringstream ss(lineContent); int wordCount = 0; - double r, theta, z; + double r, theta, z; int i,j,k; double Erho,Etheta; double Ex, Ey, Ez; @@ -124,19 +130,20 @@ namespace locust { if (wordCount == 0) { - r = std::stod(token); + r = std::stod(token); i = (int)((r-fDim1_min)/(fDim1_max-fDim1_min)*(fnPixel1)); // var1 position - } + } else if (wordCount == 1) - { - theta = std::stod(token); - j = (int)((theta-fDim2_min)/(fDim2_max-fDim2_min)*(fnPixel2)); // var2 position - } + { + theta = std::stod(token); + j = (int)((theta-fDim2_min)/(fDim2_max-fDim2_min)*(fnPixel2)); // var2 position + } else if (wordCount == 2) - { - z = std::stod(token); - k = (int)((z-fDim3_min)/(fDim3_max-fDim3_min)*(fnPixel3)); // var3 position - } + { + z = std::stod(token); + z += fZshift; + k = (int)((z-fDim3_min)/(fDim3_max-fDim3_min)*(fnPixel3)); // var3 position + } else if (wordCount == 3) Ex = std::stod(token); // mode E field value else if (wordCount == 4) Ey = std::stod(token); // mode E field value else if (wordCount == 5) Ez = std::stod(token); // mode E field value @@ -148,10 +155,10 @@ namespace locust ++wordCount; } - if ((i==fnPixel1) or (j==fnPixel2) or (k==fnPixel3)) - { - continue; - } + if ((i==fnPixel1) or (j==fnPixel2) or (k==fnPixel3)) + { + continue; + } if ((i>=fnPixel1) or (j>=fnPixel2) or (k>=fnPixel3)) { LERROR(lmclog,"Imported mode map dimensions don't agree with those in \"" << aFilename <<".\" Double check dim[1,2,3]-max."); @@ -159,26 +166,27 @@ namespace locust } //Must convert E field from cartesian coordinates to cylindrical coordinates - if(r<1.e-10) - { - Erho = 0.; - Etheta = 0.; - } - else - { + if(r<1.e-10) + { + Erho = 0.; + Etheta = 0.; + } + else + { Erho = ((Ex * r*cos(theta)) + Ey * r*sin(theta)) / r; - Etheta = ((Ey * r*cos(theta)) - Ex * r*sin(theta)) / r; - } + Etheta = ((Ey * r*cos(theta)) - Ex * r*sin(theta)) / r; + } std::vector E_input = {Erho,Etheta,Ez}; fModeMapTE_E[i][j][k] = E_input; -// printf("read var1 is %g, var2 is %g, E is %g\n", fModeMapTE_E.back()[0], fModeMapTE_E.back()[1], fModeMapTE_E.back()[2]); + //printf("read var1 is %g, var2 is %g, E is %g\n", fModeMapTE_E.back()[0], fModeMapTE_E.back()[1], fModeMapTE_E.back()[2]); + } } modeMapFile.close(); - //Reset dimensions from import file to actual cavity dimensions in case they don't match up + //Reset dimensions from import file to actual cavity dimensions in case they don't match up MatchCavityDimensions(aParam); return true; diff --git a/Source/Fields/LMCModeMapCavity.hh b/Source/Fields/LMCModeMapCavity.hh index d6f22148..894eed5d 100644 --- a/Source/Fields/LMCModeMapCavity.hh +++ b/Source/Fields/LMCModeMapCavity.hh @@ -62,6 +62,7 @@ namespace locust int fnPixel1, fnPixel2, fnPixel3; double fDim1_min, fDim2_min, fDim3_min; double fDim1_max, fDim2_max, fDim3_max; + double fZshift; std::vector< std::vector< std::vector< std::vector< double >>>> fModeMapTE_E; }; diff --git a/Source/Generators/LMCArraySignalGenerator.cc b/Source/Generators/LMCArraySignalGenerator.cc index 7f55d128..ecd3f7ac 100644 --- a/Source/Generators/LMCArraySignalGenerator.cc +++ b/Source/Generators/LMCArraySignalGenerator.cc @@ -55,6 +55,18 @@ namespace locust ArraySignalGenerator::~ArraySignalGenerator() { + if ( fAntennaElementPositioner != NULL ) + { + delete fAntennaElementPositioner; + } + if ( fTransmitter != NULL ) + { + delete fTransmitter; + } + if ( fPowerCombiner != NULL ) + { + delete fPowerCombiner; + } } void ArraySignalGenerator::SetParameters( const scarab::param_node& aParam ) @@ -347,11 +359,12 @@ namespace locust bool ArraySignalGenerator::RecordRunParameters( Signal* aSignal ) { +#ifdef ROOT_FOUND fInterface->aRunParameter = new RunParameters(); fInterface->aRunParameter->fSamplingRateMHz = fAcquisitionRate; fInterface->aRunParameter->fDecimationFactor = aSignal->DecimationFactor(); fInterface->aRunParameter->fLOfrequency = fLO_Frequency; - +#endif return true; } diff --git a/Source/Generators/LMCCavitySignalGenerator.cc b/Source/Generators/LMCCavitySignalGenerator.cc index 4a6317ca..d5a62847 100644 --- a/Source/Generators/LMCCavitySignalGenerator.cc +++ b/Source/Generators/LMCCavitySignalGenerator.cc @@ -25,12 +25,12 @@ namespace locust Generator( aName ), fDoGenerateFunc( &CavitySignalGenerator::DoGenerateTime ), fLO_Frequency( 0. ), - fDeltaT( 0. ), + fDeltaT( 0. ), gxml_filename("blank.xml"), fphiLO(0.), fNPreEventSamples( 150000 ), fRandomPreEventSamples( false ), - fTrackDelaySeed( 0 ), + fTrackDelaySeed( 0 ), fThreadCheckTime(100), fKassNeverStarted( false ), fAliasedFrequencies( false ), @@ -46,6 +46,22 @@ namespace locust CavitySignalGenerator::~CavitySignalGenerator() { + if (fPowerCombiner != NULL) + { + delete fPowerCombiner; + } + if (fFieldCalculator != NULL) + { + delete fFieldCalculator; + } + if (fTFReceiverHandler != NULL) + { + delete fTFReceiverHandler; + } + if (fAnalyticResponseFunction != NULL) + { + delete fAnalyticResponseFunction; + } } @@ -123,7 +139,6 @@ namespace locust else if (!fInterface->fbWaveguide) // Cavity config follows { // No HFSS file is present: Cavity as damped harmonic oscillator - fAnalyticResponseFunction = new DampedHarmonicOscillator(); if ( !fAnalyticResponseFunction->Configure(tParam) || !(CrossCheckCavityConfig()) ) { @@ -312,7 +327,6 @@ namespace locust bool CavitySignalGenerator::RecordRunParameters( Signal* aSignal ) { #ifdef ROOT_FOUND - fInterface->aRunParameter = new RunParameters(); fInterface->aRunParameter->fSamplingRateMHz = fAcquisitionRate; fInterface->aRunParameter->fDecimationFactor = aSignal->DecimationFactor(); fInterface->aRunParameter->fLOfrequency = fLO_Frequency; @@ -366,6 +380,11 @@ namespace locust LPROG(lmclog,"Running some cavity cross-checks ..."); + double timeResolution = 0.; + double thresholdFactor = 0.; + double cavityFrequency = 0.; + double qExpected = 0.; + for (int mu=0; muGetDHOTimeResolution(bTE, l, m, n); - double thresholdFactor = fAnalyticResponseFunction->GetDHOThresholdFactor(bTE, l, m, n); - double cavityFrequency = fAnalyticResponseFunction->GetCavityFrequency(bTE, l, m, n); - double qExpected = fAnalyticResponseFunction->GetCavityQ(bTE, l, m, n); + timeResolution = fAnalyticResponseFunction->GetDHOTimeResolution(bTE, l, m, n); + thresholdFactor = fAnalyticResponseFunction->GetDHOThresholdFactor(bTE, l, m, n); + cavityFrequency = fAnalyticResponseFunction->GetCavityFrequency(bTE, l, m, n); + qExpected = fAnalyticResponseFunction->GetCavityQ(bTE, l, m, n); aCavityUtility.SetOutputFile(fUnitTestRootFile); if (!aCavityUtility.CheckCavityQ( bTE, l, m, n, timeResolution, thresholdFactor, cavityFrequency, qExpected )) { @@ -386,6 +405,19 @@ namespace locust } } +#ifdef ROOT_FOUND + fInterface->aRunParameter = new RunParameters(); + fInterface->aRunParameter->fSimulationType = "cavity"; + if ( cavityFrequency > 20.e9 ) + { + fInterface->aRunParameter->fSimulationSubType = "cca"; + } + else + { + fInterface->aRunParameter->fSimulationSubType = "lfa"; + } +#endif + return true; } diff --git a/Source/Generators/LMCLowPassFilterFFTGenerator.cc b/Source/Generators/LMCLowPassFilterFFTGenerator.cc index 69c7dce7..9e12d87e 100644 --- a/Source/Generators/LMCLowPassFilterFFTGenerator.cc +++ b/Source/Generators/LMCLowPassFilterFFTGenerator.cc @@ -74,16 +74,19 @@ namespace locust // Count trailing zeroes in window. int LowPassFilterFFTGenerator::GetEndingMargin( Signal* aSignal, int windowsize, int nwin, int ch ) { - int endingMargin = 0; - for (int i=0; iLongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + (nwin+1)*windowsize - i ][0]) > 0.) - { - endingMargin = i; - break; - } - } - return endingMargin; + int endingMargin = 0; + for (int i=0; iLongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + (nwin+1)*windowsize - i ][0]) > 0.) + { + endingMargin = i; + break; + } + } + } + return endingMargin; } @@ -164,8 +167,10 @@ namespace locust fftw_destroy_plan(ForwardPlan); fftw_destroy_plan(ReversePlan); - delete [] SignalComplex; - delete [] FFTComplex; + fftw_free( SignalComplex ); + fftw_free( FFTComplex ); + SignalComplex = NULL; + FFTComplex = NULL; } else { @@ -175,7 +180,6 @@ namespace locust } // nwin } // NCHANNELS - return true; } diff --git a/Source/IO/LMCRunParameters.cc b/Source/IO/LMCRunParameters.cc index 6a389aa4..e2df78b2 100644 --- a/Source/IO/LMCRunParameters.cc +++ b/Source/IO/LMCRunParameters.cc @@ -17,7 +17,32 @@ ClassImp(locust::RunParameters); namespace locust { - RunParameters::RunParameters() {} + RunParameters::RunParameters() : + fNoise( 0. ), + fLOfrequency( 0. ), + fSamplingRateMHz( 0. ), + fDecimationFactor( 0. ), + fDataType( "simulation" ), + fSimulationType( "rectangular-waveguide" ), + fSimulationSubType( "none" ), + fRunID( 0 ) + { + time_t rawtime; + struct tm * timeInfo; + time (&rawtime); + timeInfo = localtime (&rawtime); + + struct timeval tv; + gettimeofday(&tv, NULL); + + int tDay = timeInfo->tm_mday; + int tMonth = timeInfo->tm_mon + 1; + int tYear = timeInfo->tm_year - 100; + int tMicrosec = tv.tv_usec; + + fRunID = 1e12 + tDay*1e10 + tMonth*1e8 + tYear*1e6 + tMicrosec; + + } RunParameters::~RunParameters() {} } diff --git a/Source/IO/LMCRunParameters.hh b/Source/IO/LMCRunParameters.hh index 867202de..1752e09b 100644 --- a/Source/IO/LMCRunParameters.hh +++ b/Source/IO/LMCRunParameters.hh @@ -15,6 +15,8 @@ #include "TObject.h" #include "LMCRunParameters.hh" +#include "time.h" +#include namespace locust { @@ -30,6 +32,10 @@ namespace locust double fLOfrequency; double fSamplingRateMHz; double fDecimationFactor; + std::string fDataType; + std::string fSimulationType; + std::string fSimulationSubType; + long int fRunID; ClassDef(RunParameters,1) // Root syntax. diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc index 3fc6db83..473b3737 100644 --- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc +++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc @@ -37,6 +37,10 @@ namespace locust } CyclotronRadiationExtractor::~CyclotronRadiationExtractor() { + if (fFieldCalculator != NULL) + { + delete fFieldCalculator; + } } bool CyclotronRadiationExtractor::Configure() @@ -62,11 +66,11 @@ namespace locust bool CyclotronRadiationExtractor::UpdateTrackProperties( Kassiopeia::KSParticle &aFinalParticle, unsigned index, bool bStart ) { +#ifdef ROOT_FOUND double tLOfrequency = fInterface->aRunParameter->fLOfrequency; // Hz double tSamplingRate = fInterface->aRunParameter->fSamplingRateMHz; // MHz double tOffset = -tLOfrequency + tSamplingRate * 1.e6 / 2.; // Hz double tTime = index / tSamplingRate / 1.e6 / fInterface->aRunParameter->fDecimationFactor; -#ifdef ROOT_FOUND if (bStart) { fStartingIndex = index; diff --git a/Source/Kassiopeia/LMCEventHold.cc b/Source/Kassiopeia/LMCEventHold.cc index bf343299..00de7638 100644 --- a/Source/Kassiopeia/LMCEventHold.cc +++ b/Source/Kassiopeia/LMCEventHold.cc @@ -8,6 +8,8 @@ #include "logger.hh" #include "LMCEventHold.hh" #include +#include + namespace locust { @@ -110,14 +112,25 @@ namespace locust { std::string tOutputPath = TOSTRING(PB_OUTPUT_DIR); std::string sFileName = tOutputPath+"/"+fTruthOutputFilename; + fJsonFileName = sFileName; + const std::string ext(".root"); + if ( sFileName != ext && + sFileName.size() > ext.size() && + sFileName.substr(sFileName.size() - ext.size()) == ".root" ) + { + // Remove .root and replace it with .json: + fJsonFileName = sFileName.substr(0, sFileName.size() - ext.size()) + ".json"; + } + else + { + LERROR(lmclog,"The output file " << fTruthOutputFilename <<"doesn't end in .root"); + exit(-1); + } #ifdef ROOT_FOUND FileWriter* aRootTreeWriter = RootTreeWriter::get_instance(); aRootTreeWriter->SetFilename(sFileName); if (fAccumulateTruthInfo) { - // TO-DO: This option should be used when running pileup. We will need to - // figure out how to explicitly increment the event structure ID, given that the - // same (identical) simulation is being run multiple times in this case. aRootTreeWriter->OpenFile("UPDATE"); } else @@ -126,11 +139,103 @@ namespace locust } aRootTreeWriter->CloseFile(); #endif + + // Open the json file: + if (!fAccumulateTruthInfo) + { + std::ofstream ost {fJsonFileName, std::ios_base::out}; + } + return true; } + bool EventHold::WriteJsonFile() + { + std::ifstream jsonFile(fJsonFileName); // Open json file for inspection + bool bNewRun = true; + std::vector v; + if (jsonFile.is_open()) + { + std::string line; + while (std::getline(jsonFile, line)) + { + LPROG( lmclog, line ); + bNewRun = !line.find("run-id"); + if (line != "}") // Avoid saving the last "}". It will be appended below. + { + v.push_back(line); + } + } + jsonFile.close(); + } + else + { + LPROG( lmclog, "A json file for meta-data was not found. Opening a new one now." ); + } + + std::ofstream ost {fJsonFileName, std::ios_base::out}; + + +#ifdef ROOT_FOUND + if (bNewRun) // If there are no run parameters in the json file yet, write them now: + { + ost << "{\n"; + ost << " \"run-id\": "<< "\"" << fInterface->aRunParameter->fRunID << "\",\n"; + ost << " \"run-parameters\": {\n"; + ost << " \"run-type\": "<< "\"" << fInterface->aRunParameter->fDataType << "\",\n"; + ost << " \"simulation-type\": "<< "\"" << fInterface->aRunParameter->fSimulationType << "\",\n"; + ost << " \"simulation-subtype\": "<< "\"" << fInterface->aRunParameter->fSimulationSubType << "\",\n"; + ost << " \"sampling-freq-mega-hz\": "<< "\"" << fInterface->aRunParameter->fSamplingRateMHz << "\"\n"; + ost << " },\n"; + } + else // otherwise re-write the file: + { + for (int i = 0; i < v.size(); i++) + { + if (i < v.size()-1) + { + ost << v[i] << "\n"; + } + else + { + ost << v[i] << ",\n"; + } + } + } + + + // Write the latest event information here: + + ost << " \"" << fInterface->anEvent->fEventID << "\": {\n"; + ost << " \"ntracks\": "<< "\"" << fInterface->anEvent->fNTracks << "\",\n"; + for (int i=0; ianEvent->fNTracks; i++) + { + ost << " \"" << fInterface->anEvent->fTrackIDs[i] << "\":\n"; + ost << " {\n"; + ost << " \"start-time\": "<< "\"" << fInterface->anEvent->fStartTimes[i] << "\",\n"; + ost << " \"end-time\": "<< "\"" << fInterface->anEvent->fEndTimes[i] << "\",\n"; + ost << " \"output-avg-frequency\": "<< "\"" << fInterface->anEvent->fOutputAvgFrequencies[i] << "\",\n"; + ost << " \"pitch-angle\": "<< "\"" << fInterface->anEvent->fPitchAngles[i] << "\",\n"; + ost << " \"avg-axial-frequency\": "<< "\"" << fInterface->anEvent->fAvgAxialFrequencies[i] << "\"\n"; + if (i < fInterface->anEvent->fNTracks-1) + { + ost << " },\n"; + } + else + { + ost << " }\n"; + } + } + ost << " }\n"; + ost << "}\n"; + +#endif + + return true; + } + bool EventHold::WriteEvent() { std::string tOutputPath = TOSTRING(PB_OUTPUT_DIR); @@ -145,6 +250,7 @@ namespace locust aRootTreeWriter->WriteRunParameters(fInterface->aRunParameter); aRootTreeWriter->CloseFile(); #endif + WriteJsonFile(); return true; } @@ -185,8 +291,11 @@ namespace locust fInterface->fEventInProgress = false; fInterface->fDigitizerCondition.notify_one(); // unlock LPROG( lmclog, "Kass is waking after event" ); +#ifdef ROOT_FOUND delete fInterface->anEvent; delete fInterface->aTrack; + delete fInterface->aRunParameter; +#endif return true; } diff --git a/Source/Kassiopeia/LMCEventHold.hh b/Source/Kassiopeia/LMCEventHold.hh index 5d490394..46bfa35d 100644 --- a/Source/Kassiopeia/LMCEventHold.hh +++ b/Source/Kassiopeia/LMCEventHold.hh @@ -31,12 +31,14 @@ namespace locust bool OpenEvent(); bool OpenFile(); bool WriteEvent(); + bool WriteJsonFile(); virtual ~EventHold(); EventHold* Clone() const; bool ConfigureByInterface(); bool Configure( const scarab::param_node& aParam ); std::string fTruthOutputFilename; + std::string fJsonFileName; bool fAccumulateTruthInfo; diff --git a/Source/Kassiopeia/LMCFieldCalculator.cc b/Source/Kassiopeia/LMCFieldCalculator.cc index cbd682cc..374cf7a6 100644 --- a/Source/Kassiopeia/LMCFieldCalculator.cc +++ b/Source/Kassiopeia/LMCFieldCalculator.cc @@ -35,6 +35,14 @@ namespace locust } FieldCalculator::~FieldCalculator() { + if (fTFReceiverHandler != NULL) + { + delete fTFReceiverHandler; + } + if (fAnalyticResponseFunction != NULL) + { + delete fAnalyticResponseFunction; + } } bool FieldCalculator::ConfigureByInterface() diff --git a/Source/Kassiopeia/LMCTrackHold.cc b/Source/Kassiopeia/LMCTrackHold.cc index 52f39c72..edd5ddf8 100644 --- a/Source/Kassiopeia/LMCTrackHold.cc +++ b/Source/Kassiopeia/LMCTrackHold.cc @@ -63,7 +63,9 @@ namespace locust bool TrackHold::ExecutePreTrackModification(Kassiopeia::KSTrack &aTrack) { +#ifdef ROOT_FOUND fInterface->aTrack->Initialize(); +#endif fInterface->fNewTrackStarting = true; double tTime = aTrack.GetInitialParticle().GetTime(); @@ -76,12 +78,15 @@ namespace locust { if ( aTrack.GetTotalSteps() > 0) { +#ifdef ROOT_FOUND + fInterface->aTrack->TrackID = fTrackCounter; fInterface->anEvent->AddTrack( fInterface->aTrack ); +#endif + fTrackCounter += 1; } double tTime = aTrack.GetFinalParticle().GetTime(); LWARN(lmclog,"LMCTrack " << fTrackCounter << " is complete at Kass time " << tTime << " with total steps " << aTrack.GetTotalSteps() ); - fTrackCounter += 1; return true; } diff --git a/Source/RxComponents/LMCDampedHarmonicOscillator.cc b/Source/RxComponents/LMCDampedHarmonicOscillator.cc index f4957aa1..5192f98c 100644 --- a/Source/RxComponents/LMCDampedHarmonicOscillator.cc +++ b/Source/RxComponents/LMCDampedHarmonicOscillator.cc @@ -19,7 +19,9 @@ namespace locust fCavityFrequencyDefault( 1.067e9 ), fCavityQDefault( 1000 ), fThresholdFactorDefault ( 0.25 ), - fHannekePowerFactorDefault( 1. ) + fHannekePowerFactorDefault( 1. ), + fNModes( 2 ), + fTFReceiverHandler( 0 ) {} DampedHarmonicOscillator::~DampedHarmonicOscillator() {} @@ -277,7 +279,7 @@ namespace locust } } // i - delete aSignal; + aSignal->Reset(); return convolutionMag; @@ -369,6 +371,7 @@ namespace locust } SetGFarray( tGFArray); // now normalized. + delete fTFReceiverHandler; if ( tGFArray.size() < 1 ) return false; else return true; diff --git a/Source/Utilities/LMCCavityUtility.cc b/Source/Utilities/LMCCavityUtility.cc index a202c2a9..23aa7ede 100644 --- a/Source/Utilities/LMCCavityUtility.cc +++ b/Source/Utilities/LMCCavityUtility.cc @@ -193,9 +193,11 @@ namespace locust #ifdef ROOT_FOUND if (fWriteOutputFile) WriteRootHisto(nSteps, freqArray, gainArray); #endif - delete aSignal; - delete freqArray; - delete gainArray; + aSignal->Reset(); + delete[] freqArray; + delete[] gainArray; + delete fTFReceiverHandler; + delete fAnalyticResponseFunction; LPROG( testlog, "\nSummary:"); LPROG( testlog, "dho-threshold-factor is " << dhoThresholdFactor );