diff --git a/offline/packages/PHGarfield/MergeGasFiles.cc b/offline/packages/PHGarfield/MergeGasFiles.cc index 7bbbc3d3bd..8007c5060d 100644 --- a/offline/packages/PHGarfield/MergeGasFiles.cc +++ b/offline/packages/PHGarfield/MergeGasFiles.cc @@ -1,222 +1,212 @@ -#include - -#include -#include -#include - #include +#include -#include #include #include -#include // for pi +#include +#include +#include #include +#include +#include +#include +#include -int main() -{ - // This is a utility to test whether the process of "merging" files is actually different from a single file. - // It may be of no further use afterthe development was complete. - // TKH 6/2/2026 - int nValid = 10000; - TNtuple *Validity = new TNtuple("Validity", "Validity", "Valid:e:b:a:Vxerr:Vyerr:Vzerr"); - - // New version chooses to NOT write output to a file (which seems broken), - // but to instead just tries to merge the files and validate the copy in memory. - const std::string dir = "gasfiles"; - - auto filename = [&](const int i) - { - return dir + "/PART_" + std::to_string(i) + ".gas"; - }; - Garfield::MediumMagboltz gas; - Garfield::MediumMagboltz gas0; +namespace fs = std::filesystem; +std::string mergedName(const std::string& path, unsigned int eindex); - const std::string first = filename(0); - if (!std::filesystem::exists(first)) - { - std::cerr << "Missing first gas file: " << first << std::endl; - return 1; - } +bool searchAndUnpackDirectory( + const std::string& directoryPath, + std::set& Eindices, + std::set& Bindices, + std::map, std::string>& FileList); - if (!gas.LoadGasFile(first)) - { - std::cerr << "Failed to load " << first << std::endl; - return 1; - } - - // Gas 0 only loads the FIRST file. This will test the memory validity of the merge... - if (!gas0.LoadGasFile(first)) - { - std::cerr << "Failed to load " << first << std::endl; - return 1; - } - - for (int i = 1;; ++i) - { - const std::string file = filename(i); - - if (!std::filesystem::exists(file)) +int main(int argc, char* argv[]) +{ + try { - std::cout << "Stopping at first missing file: " << file << std::endl; - break; + if (argc != 3) + { + std::cerr << "Usage:\n" + << argv[0] + << " path_to_gasfiles name_of_output_file\n"; + return 1; + } + + const std::string path = argv[1]; + const std::string output = path + "/" + argv[2]; + + std::set Eindices; + std::set Bindices; + std::map, std::string> FileList; + + if (!searchAndUnpackDirectory(path, Eindices, Bindices, FileList)) + { + std::cerr << PHWHERE << " Imperfect directory." << std::endl; + return 1; + } + + Garfield::MediumMagboltz gas; + + for (const auto Eindex : Eindices) + { + bool firstB = true; + + for (const auto Bindex : Bindices) + { + const auto it = FileList.find({Eindex, Bindex}); + if (it == FileList.end()) + { + std::cerr << PHWHERE << " Missing file for E=" << Eindex + << " B=" << Bindex << std::endl; + return 1; + } + + const std::string& nextfile = it->second; + + if (firstB) + { + gas.LoadGasFile(nextfile); + firstB = false; + } + else + { + gas.MergeGasFile(nextfile, true); + } + } + + const std::string mergedFile = mergedName(path, Eindex); + + std::cout << "Writing " << mergedFile << std::endl; + gas.WriteGasFile(mergedFile); + + std::vector nE; + std::vector nB; + std::vector nA; + gas.GetFieldGrid(nE, nB, nA); + + std::cout << "Merged Gas File created: "<< mergedFile + << " with Grid Dimensions: " + << nE.size() << " E-fields, " + << nB.size() << " B-fields, " + << nA.size() << " Angles." << std::endl; + } + + bool firstE = true; + + for (const auto Eindex : Eindices) + { + //if (Eindex > 10) {break;} + const std::string mergedFile = mergedName(path, Eindex); + + if (!fs::exists(mergedFile)) + { + std::cerr << PHWHERE << " Missing merged file " << mergedFile << std::endl; + return 1; + } + + if (firstE) + { + gas.LoadGasFile(mergedFile); + firstE = false; + } + else + { + gas.MergeGasFile(mergedFile, true); + } + } + + std::cout << "Writing final file " << output << std::endl; + gas.WriteGasFile(output); + + std::vector nE; + std::vector nB; + std::vector nA; + gas.GetFieldGrid(nE, nB, nA); + + std::cout << "Final Gas File created: "<< output + << " with Grid Dimensions: " + << nE.size() << " E-fields, " + << nB.size() << " B-fields, " + << nA.size() << " Angles." << std::endl; + + return 0; } - std::cout << "Merging " << file << std::endl; - - if (!gas.MergeGasFile(file, true)) + catch (const std::exception& e) { - std::cerr << "Failed to merge " << file << std::endl; + std::cerr << PHWHERE << " Exception: " << e.what() << std::endl; return 1; } - } - - // Don't write out since it crashes? - // gas.WriteGasFile("test.gas"); - - // Now perform the validation test... - double emin = 400; - double emax = 400; - // double ne=1; - - double bmin = 1.15; - double bmax = 1.45; - double nb = 50; - - double amin = 0.0; - double amax = 0.2; - // double na=50; - - // Initialize using the current system time - TRandom3 Randy; - Randy.SetSeed(PHRandomSeed()); // new initialization each run - std::cout << std::endl - << std::endl - << "Valid Calls: " << std::endl; - for (int i = 0; i < nValid; i++) - { - double eMag = Randy.Uniform(emin, emax); - double bMag = Randy.Uniform(bmin, bmin + 0.2 * (bmax - bmin) / nb); // Comes from file0... - double a = Randy.Uniform(amin, amax); - double PHI = Randy.Uniform(0.0, 2.0 * std::numbers::pi); - - double ex = 0; - double ey = 0; - double ez = eMag; - - double bx = bMag * sin(a) * cos(PHI); - double by = bMag * sin(a) * sin(PHI); - double bz = bMag * cos(a); - - double vx; - double vy; - double vz; - - double vx0; - double vy0; - double vz0; - - gas.ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz); - gas0.ElectronVelocity(ex, ey, ez, bx, by, bz, vx0, vy0, vz0); - - /* - std::cout << " i:" <Fill(1, sqrt(ex * ex + ey * ey + ez * ez), sqrt(bx * bx + by * by + bz * bz), a, DelVx, DelVy, DelVz); - } - - std::cout << std::endl - << std::endl - << "Invalid Calls: " << std::endl; - for (int i = 0; i < nValid; i++) - { - double eMag = Randy.Uniform(emin, emax); - double bMag = Randy.Uniform(bmin + 5.0 * (bmax - bmin) / nb, bmax); // Comes from beyond file0... - double a = Randy.Uniform(amin, amax); - double PHI = Randy.Uniform(0.0, 2.0 * std::numbers::pi); - - double ex = 0; - double ey = 0; - double ez = eMag; - - double bx = bMag * sin(a) * cos(PHI); - double by = bMag * sin(a) * sin(PHI); - double bz = bMag * cos(a); - - double vx; - double vy; - double vz; - double vx0; - double vy0; - double vz0; - gas.ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz); - gas0.ElectronVelocity(ex, ey, ez, bx, by, bz, vx0, vy0, vz0); - /* - std::cout << " i:" <Fill(0, sqrt(ex * ex + ey * ey + ez * ez), sqrt(bx * bx + by * by + bz * bz), a, DelVx, DelVy, DelVz); - } + catch (...) + { + std::cerr << PHWHERE << " Unknown exception." << std::endl; + return 1; + } +} - TFile *output = new TFile("GarfieldValidity.root", "RECREATE"); - Validity->Write(); - output->Close(); +bool searchAndUnpackDirectory(const std::string& directoryPath, std::set &Eindices, std::set &Bindices, std::map, std::string> &FileList) +{ + // Check if the directory exists and is valid + if (!fs::exists(directoryPath) || !fs::is_directory(directoryPath)) + { + std::cerr << "Error: Invalid directory path." << std::endl; + return false; + } + + std::regex filePattern(R"(^E([0-9]{3})_B([0-9]{3})\.gas$)"); + std::smatch matchResults; + + // Iterate through all items in the directory + for (const auto& entry : fs::directory_iterator(directoryPath)) + { + // Only process regular files + if (entry.is_regular_file()) + { + std::string filename = entry.path().filename().string(); + + // Check if the filename matches our target pattern + if (std::regex_match(filename, matchResults, filePattern)) + { + // matchResults[1] contains the string after 'E' + // matchResults[2] contains the string after 'B' + // std::stoul automatically handles leading zeros + unsigned int eValue = std::stoul(matchResults[1].str()); + unsigned int bValue = std::stoul(matchResults[2].str()); + Eindices.insert(eValue); + Bindices.insert(bValue); + FileList[{eValue, bValue}] = entry.path().string(); + } + } + } + + // Validate the results. + if (Eindices.empty()) { return false; } + if (Bindices.empty()) { return false; } + + unsigned int maxE = *Eindices.rbegin(); + unsigned int maxB = *Bindices.rbegin(); + for (unsigned int i=0; i<=maxE; i++) + { + for (unsigned int j=0; j<=maxB; j++) + { + if ( !FileList.contains({i,j}) ) { return false; } + } + } + + std::cout << " *** Gas File List Valid ***" << std::endl; + std::cout << "Electric field indices 0 --> " << maxE << std::endl; + std::cout << "Magnetic field indices 0 --> " << maxB << std::endl; + + return true; +} - return 0; +std::string mergedName(const std::string& path, unsigned int eindex) +{ + std::ostringstream name; + name << path << "/MERGED_E" + << std::setw(3) << std::setfill('0') << eindex + << ".gas"; + return name.str(); } diff --git a/offline/packages/PHGarfield/PHGarfield.cc b/offline/packages/PHGarfield/PHGarfield.cc index 157ffb893d..259c8803dc 100644 --- a/offline/packages/PHGarfield/PHGarfield.cc +++ b/offline/packages/PHGarfield/PHGarfield.cc @@ -1,5 +1,5 @@ #include "PHGarfield.h" - +#include #include #include @@ -20,12 +20,30 @@ #include #include // for basic_ostream, operat... #include +#include +#include +#include +#include +#include + +namespace fs = std::filesystem; PHGarfield::PHGarfield(const std::string& name) - : SubsysReco(name) + : SubsysReco(name), + //m_defaultGasfile("/sphenix/user/hemmick/gasfiles_20260624/Ar75_CF20_iso5.gas") + m_defaultGasfile("/sphenix/user/hemmick/gasfiles_20260624") { } +PHGarfield::~PHGarfield() +{ + // Housekeeping. + delete m_field; + delete m_cdbTPCMAPttree; + delete m_component; + delete m_gas; +} + int PHGarfield::InitRun(PHCompositeNode* /*topNode*/) { if (Verbosity() > 1) @@ -49,7 +67,16 @@ int PHGarfield::InitRun(PHCompositeNode* /*topNode*/) { GetMagneticFieldTesla(x, y, z, bx, by, bz); }); m_component->SetElectricField([this](double x, double y, double z, double& ex, double& ey, double& ez) { GetElectricFieldVcm(x, y, z, ex, ey, ez); }); - InitializeGas("/direct/phenix+u/workarea/hemmick/code.sphenix/tkh/gas/gasfiles/"); + + // Here we fetch the gas from the CDB + std::string gasfile = m_cdb->getUrl("PHGARFIELD_GAS"); + if (gasfile.empty() || !fs::exists(gasfile)) + { + std::cerr << PHWHERE << " Missing CDB gasfile: " << gasfile << std::endl; + std::cerr << PHWHERE << " Using default gasfile: " << m_defaultGasfile << std::endl; + gasfile = m_defaultGasfile; + } + InitializeGas(gasfile); // Diagnostic during code development... FillRadii(); @@ -113,6 +140,25 @@ void PHGarfield::PrintGarfield(double x, double y, double z) const << std::endl; } +void PHGarfield::PrintGasSummary() const +{ + if (!m_GasFilesLoaded) + { + std::cerr << PHWHERE << "No Gas File(s) have been successfully loaded." << std::endl; + return; + } + + std::vector nE; + std::vector nB; + std::vector nA; + m_gas->GetFieldGrid(nE, nB, nA); + + std::cout << "Gas File Grid Dimensions: " << std::endl; + std::cout << nE.size() << " E-fields ranging from " << nE.front() << " to " << nE.back() << std::endl; + std::cout << nB.size() << " B-fields ranging from " << nB.front() << " to " << nB.back() << std::endl; + std::cout << nA.size() << " Angles ranging from " << nA.front() << " to " << nA.back() << std::endl; +} + void PHGarfield::PrintMaps() const { // Print out a few test points of the Garfield information @@ -195,52 +241,77 @@ void PHGarfield::GetElectricFieldVcm(double x_cm, double y_cm, double z_cm, doub ez_vcm = z_cm > 0 ? -400.0 : 400.0; } -void PHGarfield::InitializeGas(const std::string &dir) +void PHGarfield::InitializeGas(const std::string &name) { // Create and fill the gas object so that we can trace particles through the gas... m_gas = new Garfield::MediumMagboltz(); - auto filename = [&](const int i) - { return dir + "/PART_" + std::to_string(i) + ".gas"; }; - - const std::string first = filename(0); - if (!std::filesystem::exists(first)) + if (!std::filesystem::exists(name)) { - std::cerr << "Missing first gas file: " << first << std::endl; + std::cerr << "Missing gas file or gas directory: " << name << std::endl; return; } - if (!m_gas->LoadGasFile(first)) - { - std::cerr << "Failed to load " << first << std::endl; - return; - } - - for (int i = 1;; ++i) - { - const std::string file = filename(i); - - if (!std::filesystem::exists(file)) + if (fs::is_regular_file(name)) { - std::cout << "Stopping at first missing file: " << file << std::endl; - break; + std::cout << "Loading Garfield gas from file: " << name << std::endl; + if (!m_gas->LoadGasFile(name)) + { + std::cerr << "Failed to load " << name << std::endl; + return; + } + m_GasFilesLoaded = true; } - - std::cout << "Merging " << file << std::endl; - - if (!m_gas->MergeGasFile(file, true)) + else if (fs::is_directory(name)) { - std::cerr << "Failed to merge " << file << std::endl; - return; + std::cout << "Loading Garfield gas from directory: " << name << std::endl; + std::regex filePattern(R"(^MERGED_E([0-9]{3})\.gas$)"); + std::smatch matchResults; + + // Iterate through all items in the directory + // NOTE: Map assures that files are properly ordered when merged... + std::map FilesToMerge; + for (const auto& entry : fs::directory_iterator(name)) + { + // Only process regular files + if (entry.is_regular_file()) + { + std::string filepath = entry.path().string(); + std::string filename = entry.path().filename().string(); + + // Check if the filename matches our target pattern + if (std::regex_match(filename, matchResults, filePattern)) + { + //std::cout << "matchResults: " << matchResults[1].str() << std::endl; + FilesToMerge[std::stoul( matchResults[1].str() )]=filepath; + } + } + } + bool firstE = true; + for (const auto& [key, filepath] : FilesToMerge) + { + if (firstE) + { + m_gas->LoadGasFile(filepath); + firstE = false; + m_GasFilesLoaded = true; + } + else + { + m_gas->MergeGasFile(filepath, true); + m_GasFilesLoaded = true; + } + } } - } + + PrintGasSummary(); } -int PHGarfield::process_event(PHCompositeNode*) +int PHGarfield::process_event(PHCompositeNode* topNode) { // Initial implementation doesn't do anything event-by-event. // Nonetheless, a future user might want do do something here... - + (void) topNode; return Fun4AllReturnCodes::EVENT_OK; } diff --git a/offline/packages/PHGarfield/PHGarfield.h b/offline/packages/PHGarfield/PHGarfield.h index 47d3118898..ff2d8bb957 100644 --- a/offline/packages/PHGarfield/PHGarfield.h +++ b/offline/packages/PHGarfield/PHGarfield.h @@ -21,15 +21,16 @@ class PHGarfield : public SubsysReco { public: PHGarfield(const std::string &name = "PHGarfield"); - ~PHGarfield() override = default; + ~PHGarfield() override; int InitRun(PHCompositeNode *) override; - int process_event(PHCompositeNode * /*topNode*/) override; + int process_event(PHCompositeNode * topNode) override; bool StopHere(const double x, const double y, const double z, const double zPrevious); void PrintMaps() const; void PrintGarfield(double x, double y, double z) const; + void PrintGasSummary() const; // These are left in public namespace for easy plotting macros... // The user is encouraged to add more routine to fit their analysis goals... @@ -40,14 +41,16 @@ class PHGarfield : public SubsysReco private: void GetMagneticFieldTesla(double x_cm, double y_cm, double z_cm, double &bx_t, double &by_t, double &bz_t) const; // Feeds magnetic field to Garfield void GetElectricFieldVcm(double x_cm, double y_cm, double z_cm, double &ex_vcm, double &ey_vcm, double &ez_vcm) const; // Feeds electric field to Garfield - void InitializeGas(const std::string &dir); + void InitializeGas(const std::string &name); // Acepts a file or a directory void FillRadii(); static double bounder(double phi, double phi_min); CDBTTree *m_cdbTPCMAPttree{nullptr}; // Locations of the pads from CDB... - PHField3DCartesian *m_field{nullptr}; // The stanards sPHENIX field holding container. + PHField3DCartesian *m_field{nullptr}; // The standard sPHENIX field holding container. Garfield::ComponentUser *m_component{nullptr}; // This handles the interface of the electric and magnetic fields as handed to Garfield Garfield::MediumMagboltz *m_gas{nullptr}; // This is the pre-tabulated gas properties required by Garfield... + std::string m_defaultGasfile; + bool m_GasFilesLoaded{false}; // These are utilities for a spot check of the overall routine: // std::string calibdir;