packages/PHGarfield/PHGarfield.cc

0001 #include "PHGarfield.h"
0002 #include <phool/phool.h>
0003 #include <cdbobjects/CDBTTree.h>
0004
0005 #include <phfield/PHField3DCartesian.h>
0006
0007 #include <ffamodules/CDBInterface.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010
0011 #include <TPolyLine3D.h>
0012
0013 #include <CLHEP/Units/SystemOfUnits.h>
0014
0015 #include <Garfield/ComponentUser.hh>
0016 #include <Garfield/MediumMagboltz.hh>
0017
0018 #include <cmath>
0019 #include <filesystem>
0020 #include <functional>
0021 #include <iostream>  // for basic_ostream, operat...
0022 #include <vector>
0023 #include <regex>
0024 #include <string>
0025 #include <utility>
0026 #include <sstream>
0027 #include <iomanip>
0028
0029 namespace fs = std::filesystem;
0030
0031 PHGarfield::PHGarfield(const std::string& name)
0032   : SubsysReco(name),
0033     //m_defaultGasfile("/sphenix/user/hemmick/gasfiles_20260624/Ar75_CF20_iso5.gas")
0034     m_defaultGasfile("/sphenix/user/hemmick/gasfiles_20260624")
0035 {
0036 }
0037
0038 PHGarfield::~PHGarfield()
0039 {
0040   //  Housekeeping.
0041   delete m_field;
0042   delete m_cdbTPCMAPttree;
0043   delete m_component;
0044   delete m_gas;
0045 }
0046
0047 int PHGarfield::InitRun(PHCompositeNode* /*topNode*/)
0048 {
0049   if (Verbosity() > 1)
0050   {
0051     std::cout << "PHGarfield::InitRun(PHCompositeNode *topNode) Initializing" << std::endl;
0052   }
0053   CDBInterface* m_cdb = CDBInterface::instance();
0054
0055   //  Here we use the CDBInterface to set up the magnetic field map:
0056   std::string url = m_cdb->getUrl("FIELDMAP_TRACKING");
0057   m_field = new PHField3DCartesian(url, 1.0);
0058
0059   //  Here we use the CDBInterface to set up the channel making of the TPC:
0060   std::string text = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP");
0061   m_cdbTPCMAPttree = new CDBTTree(text);
0062   m_cdbTPCMAPttree->LoadCalibrations();
0063
0064   //  Make the Garfield Component and register the methods that will interface to our fields...
0065   m_component = new Garfield::ComponentUser();
0066   m_component->SetMagneticField([this](double x, double y, double z, double& bx, double& by, double& bz)
0067                                 { GetMagneticFieldTesla(x, y, z, bx, by, bz); });
0068   m_component->SetElectricField([this](double x, double y, double z, double& ex, double& ey, double& ez)
0069                                 { GetElectricFieldVcm(x, y, z, ex, ey, ez); });
0070
0071   // Here we fetch the gas from the CDB
0072   std::string gasfile = m_cdb->getUrl("PHGARFIELD_GAS");
0073   if (gasfile.empty() || !fs::exists(gasfile))
0074     {
0075       std::cerr << PHWHERE << " Missing CDB gasfile: " << gasfile << std::endl;
0076       std::cerr << PHWHERE << " Using default gasfile: " << m_defaultGasfile << std::endl;
0077       gasfile = m_defaultGasfile;
0078     }
0079   InitializeGas(gasfile);
0080
0081   //  Diagnostic during code development...
0082   FillRadii();
0083   if (Verbosity() > 1)
0084   {
0085     PrintMaps();
0086   }
0087   return Fun4AllReturnCodes::EVENT_OK;
0088 }
0089
0090 void PHGarfield::FillRadii()
0091 {
0092   //  Unload the pad map to get the radii in a handy location:
0093   for (unsigned int side = 0; side < 2; side++)
0094   {
0095     for (unsigned int sector = 0; sector < 12; sector++)
0096     {
0097       for (unsigned int fee = 0; fee < 26; fee++)
0098       {
0099         for (unsigned int channel = 0; channel < 256; channel++)
0100         {
0101           unsigned int key = (256 * (fee)) + channel;
0102           int layer = m_cdbTPCMAPttree->GetIntValue(key, "layer");
0103           double r = m_cdbTPCMAPttree->GetDoubleValue(key, "R") / CLHEP::cm;
0104           if (layer > 6)
0105           {
0106             radii[layer - 7] = r;
0107           }
0108         }
0109       }
0110     }
0111   }
0112 }
0113
0114 void PHGarfield::PrintGarfield(double x, double y, double z) const
0115 {
0116   double ex;
0117   double ey;
0118   double ez;
0119   double bx;
0120   double by;
0121   double bz;
0122   double vx;
0123   double vy;
0124   double vz;
0125   GetElectricFieldVcm(x, y, z, ex, ey, ez);
0126   GetMagneticFieldTesla(x, y, z, bx, by, bz);
0127   m_gas->ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz);
0128   std::cout << " x:" << x
0129             << " y:" << y
0130             << " z:" << z
0131             << " ex:" << ex
0132             << " ey:" << ey
0133             << " ez:" << ez
0134             << " bx:" << bx
0135             << " by:" << by
0136             << " bz:" << bz
0137             << " vx:" << vx
0138             << " vy:" << vy
0139             << " vz:" << vz
0140             << std::endl;
0141 }
0142
0143 void PHGarfield::PrintGasSummary() const
0144 {
0145   if (!m_GasFilesLoaded)
0146     {
0147       std::cerr << PHWHERE << "No Gas File(s) have been successfully loaded." << std::endl;
0148       return;
0149     }
0150
0151   std::vector<double> nE;
0152   std::vector<double> nB;
0153   std::vector<double> nA;
0154   m_gas->GetFieldGrid(nE, nB, nA);
0155
0156   std::cout << "Gas File Grid Dimensions: " << std::endl;
0157   std::cout << nE.size() << " E-fields ranging from " << nE.front() << " to " << nE.back() << std::endl;
0158   std::cout << nB.size() << " B-fields ranging from " << nB.front() << " to " << nB.back() << std::endl;
0159   std::cout << nA.size() << " Angles   ranging from " << nA.front() << " to " << nA.back() << std::endl;
0160 }
0161
0162 void PHGarfield::PrintMaps() const
0163 {
0164   //  Print out a few test points of the Garfield information
0165   PrintGarfield(0.0, 0.0, 0.1);
0166   PrintGarfield(0.0, 0.0, 100.0);
0167   PrintGarfield(0.0, 40.0, 100.1);
0168   PrintGarfield(0.0, 78.0, 010.1);
0169
0170   //  Print out the pad coordinate map:
0171   int MAX = 10;
0172   int prints = 0;
0173   for (unsigned int side = 0; side < 2; side++)
0174   {
0175     for (unsigned int sector = 0; sector < 12; sector++)
0176     {
0177       for (unsigned int fee = 0; fee < 26; fee++)
0178       {
0179         for (unsigned int channel = 0; channel < 256; channel++)
0180         {
0181           unsigned int key = (256 * (fee)) + channel;
0182           int layer = m_cdbTPCMAPttree->GetIntValue(key, "layer");
0183           double phi = ((side == 1 ? 1 : -1) * (m_cdbTPCMAPttree->GetDoubleValue(key, "phi") - std::numbers::pi / 2.)) + ((sector % 12) * std::numbers::pi / 6);
0184           double r = m_cdbTPCMAPttree->GetDoubleValue(key, "R") / CLHEP::cm;
0185
0186           phi = bounder(phi, PHI_MIN);
0187
0188           if (layer > 6)
0189           {
0190             if (prints < MAX)
0191             {
0192               prints++;
0193               std::cout << " side: " << side;
0194               std::cout << " sector: " << sector;
0195               std::cout << " fee: " << fee;
0196               std::cout << " channel: " << channel;
0197               std::cout << " layer: " << layer;
0198               std::cout << " phi: " << phi;
0199               std::cout << " r: " << r;
0200               std::cout << std::endl;
0201             }
0202           }
0203         }
0204       }
0205     }
0206   }
0207 }
0208
0209 void PHGarfield::GetMagneticFieldTesla(double x_cm, double y_cm, double z_cm, double& bx_t, double& by_t, double& bz_t) const
0210 {
0211   // NOTE:  Garfield uses  cm, V/cm, and Tesla.
0212   //        CLHEP    uses  mm, V/mm, and kiloTesla
0213   //        PHField3DCartesian follows the CLHEP conventions for magnetic fields.
0214
0215   double point[4] =
0216       {
0217           x_cm * CLHEP::cm,
0218           y_cm * CLHEP::cm,
0219           z_cm * CLHEP::cm,
0220           //(z_cm-20.0) * CLHEP::cm,
0221           0.0};
0222
0223   double bfield[3] = {0.0, 0.0, 0.0};
0224
0225   //  Get the magnetic field via the PHField3DCartesian object constructed usinf the CDB url reference.
0226   m_field->GetFieldValue(point, bfield);
0227
0228   bx_t = bfield[0] / CLHEP::tesla;
0229   by_t = bfield[1] / CLHEP::tesla;
0230   bz_t = bfield[2] / CLHEP::tesla;
0231 }
0232
0233 void PHGarfield::GetElectricFieldVcm(double x_cm, double y_cm, double z_cm, double& ex_vcm, double& ey_vcm, double& ez_vcm) const
0234 {
0235   // NOTE:  Garfield uses  cm, V/cm, and Tesla.
0236   (void) x_cm;
0237   (void) y_cm;
0238
0239   ex_vcm = 0.0;
0240   ey_vcm = 0.0;
0241   ez_vcm = z_cm > 0 ? -400.0 : 400.0;
0242 }
0243
0244 void PHGarfield::InitializeGas(const std::string &name)
0245 {
0246   //  Create and fill the gas object so that we can trace particles through the gas...
0247   m_gas = new Garfield::MediumMagboltz();
0248
0249   if (!std::filesystem::exists(name))
0250   {
0251     std::cerr << "Missing gas file or gas directory: " << name << std::endl;
0252     return;
0253   }
0254
0255   if (fs::is_regular_file(name))
0256     {
0257       std::cout << "Loading Garfield gas from file: " << name << std::endl;
0258       if (!m_gas->LoadGasFile(name))
0259     {
0260       std::cerr << "Failed to load " << name << std::endl;
0261       return;
0262     }
0263       m_GasFilesLoaded = true;
0264     }
0265   else if (fs::is_directory(name))
0266     {
0267       std::cout << "Loading Garfield gas from directory: " << name << std::endl;
0268       std::regex filePattern(R"(^MERGED_E([0-9]{3})\.gas$)");
0269       std::smatch matchResults;
0270
0271       // Iterate through all items in the directory
0272       // NOTE:  Map assures that files are properly ordered when merged...
0273       std::map<unsigned int, std::string> FilesToMerge;
0274       for (const auto& entry : fs::directory_iterator(name))
0275     {
0276       // Only process regular files
0277       if (entry.is_regular_file())
0278         {
0279           std::string filepath = entry.path().string();
0280           std::string filename = entry.path().filename().string();
0281
0282           // Check if the filename matches our target pattern
0283           if (std::regex_match(filename, matchResults, filePattern))
0284         {
0285           //std::cout << "matchResults: " << matchResults[1].str() << std::endl;
0286           FilesToMerge[std::stoul( matchResults[1].str() )]=filepath;
0287         }
0288         }
0289     }
0290       bool firstE = true;
0291       for (const auto& [key, filepath] : FilesToMerge)
0292     {
0293       if (firstE)
0294         {
0295           m_gas->LoadGasFile(filepath);
0296           firstE = false;
0297           m_GasFilesLoaded = true;
0298         }
0299       else
0300         {
0301           m_gas->MergeGasFile(filepath, true);
0302           m_GasFilesLoaded = true;
0303         }
0304     }
0305     }
0306
0307   PrintGasSummary();
0308 }
0309
0310 int PHGarfield::process_event(PHCompositeNode* topNode)
0311 {
0312   // Initial implementation doesn't do anything event-by-event.
0313   // Nonetheless, a future user might want do do something here...
0314   (void) topNode;
0315   return Fun4AllReturnCodes::EVENT_OK;
0316 }
0317
0318 double PHGarfield::bounder(double phi, double phi_min)
0319 {
0320   double phi_max = phi_min + 2.0 * std::numbers::pi;
0321   while (phi < phi_min)
0322   {
0323     phi = phi + 2.0 * std::numbers::pi;
0324   }
0325   while (phi >= phi_max)
0326   {
0327     phi = phi - 2.0 * std::numbers::pi;
0328   }
0329
0330   return phi;
0331 }
0332
0333 TPolyLine3D* PHGarfield::ReverseDrift(double x, double y, double z, double step_ns)
0334 {
0335   std::vector<double> xlist;
0336   std::vector<double> ylist;
0337   std::vector<double> zlist;
0338
0339   xlist.push_back(x);
0340   ylist.push_back(y);
0341   zlist.push_back(z);
0342
0343   double ex;
0344   double ey;
0345   double ez;
0346   double bx;
0347   double by;
0348   double bz;
0349   double vx;
0350   double vy;
0351   double vz;
0352
0353   double zPrevious = z;
0354   while (!StopHere(x, y, z, zPrevious))
0355   {
0356     zPrevious = z;
0357     GetMagneticFieldTesla(x, y, z, bx, by, bz);
0358     GetElectricFieldVcm(x, y, z, ex, ey, ez);
0359     m_gas->ElectronVelocity(ex, ey, ez, bx, by, bz, vx, vy, vz);
0360
0361     x = x - vx * step_ns;
0362     y = y - vy * step_ns;
0363     z = z - vz * step_ns;
0364
0365     xlist.push_back(x);
0366     ylist.push_back(y);
0367     zlist.push_back(z);
0368   }
0369
0370   TPolyLine3D* poly = new TPolyLine3D(xlist.size() - 1);
0371   for (unsigned int i = 0; i < xlist.size() - 1; i++)
0372   {
0373     poly->SetPoint(i, xlist[i], ylist[i], zlist[i]);
0374   }
0375
0376   return poly;
0377 }
0378
0379 bool PHGarfield::StopHere(const double x, const double y, const double z,
0380                           const double zPrevious)
0381 {
0382   const double r = std::hypot(x, y);
0383
0384   if (r < 18.0)
0385   {
0386     return true;
0387   }
0388   if (r > 82.0)
0389   {
0390     return true;
0391   }
0392   if (z > 120.0)
0393   {
0394     return true;
0395   }
0396   if (z < -120.0)
0397   {
0398     return true;
0399   }
0400
0401   // z crossed the central membrane.
0402   if (z * zPrevious < 0.0)
0403   {
0404     return true;
0405   }
0406
0407   return false;
0408 }