g4simulation/g4detectors/PHG4GeantinoIonization.cc

0001 #include "PHG4GeantinoIonization.h"
0002
0003 #include <g4main/PHG4Hit.h>
0004 #include <g4main/PHG4HitContainer.h>
0005 #include <g4main/PHG4Particle.h>
0006 #include <g4main/PHG4TruthInfoContainer.h>
0007
0008 #include <fun4all/Fun4AllReturnCodes.h>
0009
0010 #include <phparameter/PHParameters.h>
0011 #include <phparameter/PHParametersContainer.h>
0012
0013 #include <pdbcalbase/PdbParameterMapContainer.h>
0014
0015 #include <phool/PHCompositeNode.h>
0016 #include <phool/PHDataNode.h>
0017 #include <phool/PHNodeIterator.h>
0018 #include <phool/getClass.h>
0019 #include <phool/phool.h>
0020
0021 #include <cmath>
0022 #include <cstddef>
0023 #include <iostream>
0024
0025 namespace
0026 {
0027   // Pure-gas MIP stopping powers in keV/cm. They match the values used by
0028   // the current TPC and Micromegas hit-reconstruction modules.
0029   constexpr double neonMipDedx = 1.56;
0030   constexpr double argonMipDedx = 2.44;
0031   constexpr double cf4MipDedx = 7.00;
0032   constexpr double nitrogenMipDedx = 2.127;
0033   constexpr double isobutaneMipDedx = 5.93;
0034
0035   // Mean silicon MIP stopping powers in GeV/cm. The MVTX value corresponds
0036   // to 9.6 keV in 25 microns; the INTT value is the value documented by its
0037   // hit reconstruction.
0038   double mvtxMipDedx = 0.00384;
0039   double inttMipDedx = 0.00387;
0040
0041   // PHG4MicromegasDetector and PHG4MicromegasHitReco both use a fixed
0042   // Ar/isobutane 90/10 gas mixture.
0043   constexpr double tpotMipDedx =
0044       1e-6 * (0.9 * argonMipDedx + 0.1 * isobutaneMipDedx);
0045 }  // namespace
0046
0047 PHG4GeantinoIonization::PHG4GeantinoIonization(const std::string& name)
0048   : SubsysReco(name)
0049   , m_detectorConfigs{{
0050         {DetectorId::mvtx, "MVTX", "G4HIT_MVTX", true},
0051         {DetectorId::intt, "INTT", "G4HIT_INTT", true},
0052         {DetectorId::tpc, "TPC", "G4HIT_TPC", true},
0053         {DetectorId::tpot, "MICROMEGAS", "G4HIT_MICROMEGAS", true}}}
0054 {
0055 }
0056
0057 int PHG4GeantinoIonization::InitRun(PHCompositeNode* topNode)
0058 {
0059   if (!m_detectorConfigs[2].enabled)
0060   {
0061     return Fun4AllReturnCodes::EVENT_OK;
0062   }
0063
0064   // Read the gas fractions from the TPC geometry parameters, as is done in
0065   // PHG4TpcElectronDrift. This keeps the synthetic ionization consistent with
0066   // the geometry built by the macro or loaded from the CDB.
0067   auto* tpcParamsContainer =
0068       findNode::getClass<PHParametersContainer>(topNode, "G4GEO_TPC");
0069   if (!tpcParamsContainer)
0070   {
0071     // Tracking geometry loaded from the CDB initially provides the serialized
0072     // RUN-node parameters. Rebuild G4GEO_TPC from them before TPC hit
0073     // reconstruction, following PHG4TpcElectronDrift::InitRun.
0074     auto* tpcPdbParams =
0075         findNode::getClass<PdbParameterMapContainer>(topNode, "G4GEOPARAM_TPC");
0076     if (!tpcPdbParams)
0077     {
0078       std::cout << PHWHERE
0079                 << " Missing both G4GEO_TPC and G4GEOPARAM_TPC"
0080                 << std::endl;
0081       return Fun4AllReturnCodes::ABORTRUN;
0082     }
0083
0084     PHNodeIterator topIter(topNode);
0085     auto* parNode = dynamic_cast<PHCompositeNode*>(
0086         topIter.findFirst("PHCompositeNode", "PAR"));
0087     if (!parNode)
0088     {
0089       std::cout << PHWHERE << " Missing PAR node" << std::endl;
0090       return Fun4AllReturnCodes::ABORTRUN;
0091     }
0092
0093     PHNodeIterator parIter(parNode);
0094     auto* parTpcNode = dynamic_cast<PHCompositeNode*>(
0095         parIter.findFirst("PHCompositeNode", "TPC"));
0096     if (!parTpcNode)
0097     {
0098       parTpcNode = new PHCompositeNode("TPC");
0099       parNode->addNode(parTpcNode);
0100     }
0101
0102     tpcParamsContainer = new PHParametersContainer("TPC");
0103     tpcParamsContainer->CreateAndFillFrom(tpcPdbParams, "TPC");
0104     parTpcNode->addNode(
0105         new PHDataNode<PHParametersContainer>(
0106             tpcParamsContainer, "G4GEO_TPC"));
0107   }
0108
0109   const PHParameters* tpcParams = tpcParamsContainer->GetParameters(0);
0110   if (!tpcParams)
0111   {
0112     std::cout << PHWHERE << " Missing TPC geometry parameters" << std::endl;
0113     return Fun4AllReturnCodes::ABORTRUN;
0114   }
0115
0116   const double neonFraction = tpcParams->get_double_param("Ne_frac");
0117   const double argonFraction = tpcParams->get_double_param("Ar_frac");
0118   const double cf4Fraction = tpcParams->get_double_param("CF4_frac");
0119   const double nitrogenFraction = tpcParams->get_double_param("N2_frac");
0120   const double isobutaneFraction = tpcParams->get_double_param("isobutane_frac");
0121
0122   m_tpcMipDedx =
0123       1e-6 * (neonFraction * neonMipDedx +
0124               argonFraction * argonMipDedx +
0125               cf4Fraction * cf4MipDedx +
0126               nitrogenFraction * nitrogenMipDedx +
0127               isobutaneFraction * isobutaneMipDedx);
0128
0129   if (Verbosity() > 0)
0130   {
0131     std::cout << Name()
0132               << " TPC gas fractions (Ne/Ar/CF4/N2/isobutane): "
0133               << neonFraction << "/" << argonFraction << "/"
0134               << cf4Fraction << "/" << nitrogenFraction << "/"
0135               << isobutaneFraction
0136               << ", MIP dE/dx: " << m_tpcMipDedx << " GeV/cm"
0137               << std::endl;
0138   }
0139
0140   return Fun4AllReturnCodes::EVENT_OK;
0141 }
0142
0143 void PHG4GeantinoIonization::set_mvtx_mip_dedx(const double value)
0144 {
0145   mvtxMipDedx = value;
0146 }
0147
0148 void PHG4GeantinoIonization::set_intt_mip_dedx(const double value)
0149 {
0150   inttMipDedx = value;
0151 }
0152
0153 double PHG4GeantinoIonization::mip_dedx(const DetectorId detector) const
0154 {
0155   switch (detector)
0156   {
0157   case DetectorId::mvtx:
0158     return mvtxMipDedx;
0159   case DetectorId::intt:
0160     return inttMipDedx;
0161   case DetectorId::tpc:
0162     return m_tpcMipDedx;
0163   case DetectorId::tpot:
0164     return tpotMipDedx;
0165   }
0166
0167   return 0;
0168 }
0169
0170 int PHG4GeantinoIonization::process_event(PHCompositeNode* topNode)
0171 {
0172   const auto* truthInfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0173   if (!truthInfo)
0174   {
0175     std::cout << PHWHERE << " Missing G4TruthInfo" << std::endl;
0176     return Fun4AllReturnCodes::ABORTEVENT;
0177   }
0178
0179   for (const auto& config : m_detectorConfigs)
0180   {
0181     if (!config.enabled)
0182     {
0183       continue;
0184     }
0185
0186     auto* hits = findNode::getClass<PHG4HitContainer>(topNode, config.hitNodeName);
0187     if (!hits)
0188     {
0189       if (Verbosity() > 1)
0190       {
0191         std::cout << PHWHERE << " Missing optional node "
0192                   << config.hitNodeName << std::endl;
0193       }
0194       continue;
0195     }
0196
0197     process_detector(hits, truthInfo, config);
0198   }
0199
0200   return Fun4AllReturnCodes::EVENT_OK;
0201 }
0202
0203 void PHG4GeantinoIonization::process_detector(
0204     PHG4HitContainer* hits,
0205     const PHG4TruthInfoContainer* truthInfo,
0206     const DetectorConfig& config) const
0207 {
0208   std::size_t inspected = 0;
0209   std::size_t modified = 0;
0210   std::size_t missingParticle = 0;
0211   std::size_t invalidPath = 0;
0212
0213   const auto hitRange = hits->getHits();
0214   for (auto hitIter = hitRange.first; hitIter != hitRange.second; ++hitIter)
0215   {
0216     auto* hit = hitIter->second;
0217     if (!hit)
0218     {
0219       continue;
0220     }
0221
0222     ++inspected;
0223
0224     const auto* particle = truthInfo->GetParticle(hit->get_trkid());
0225     if (!particle)
0226     {
0227       ++missingParticle;
0228       continue;
0229     }
0230
0231     if (particle->get_name() != m_particleName)
0232     {
0233       continue;
0234     }
0235
0236     // Keep the operation idempotent. Current stepping actions store negative
0237     // edep/eion sentinels for geantinos. A finite nonnegative value means this
0238     // hit has already been processed.
0239     const double edep = hit->get_edep();
0240     const double eion = hit->get_eion();
0241     if (std::isfinite(edep) && edep >= 0 &&
0242         std::isfinite(eion) && eion >= 0)
0243     {
0244       continue;
0245     }
0246
0247     const double dx = hit->get_x(1) - hit->get_x(0);
0248     const double dy = hit->get_y(1) - hit->get_y(0);
0249     const double dz = hit->get_z(1) - hit->get_z(0);
0250     const double pathLength = std::sqrt(dx * dx + dy * dy + dz * dz);
0251
0252     const double mipDedx = mip_dedx(config.detector);
0253     if (!std::isfinite(pathLength) || pathLength <= 0 ||
0254         !std::isfinite(mipDedx) || mipDedx <= 0)
0255     {
0256       ++invalidPath;
0257       continue;
0258     }
0259
0260     const double syntheticEnergyDeposit = mipDedx * pathLength;
0261     const double syntheticIonization = syntheticEnergyDeposit;
0262     hit->set_edep(syntheticEnergyDeposit);
0263     hit->set_eion(syntheticIonization);
0264     ++modified;
0265
0266     if (Verbosity() > 2)
0267     {
0268       std::cout << Name() << " " << config.name
0269                 << " hit " << hitIter->first
0270                 << " track " << hit->get_trkid()
0271                 << " path length " << pathLength << " cm"
0272                 << " synthetic edep " << syntheticEnergyDeposit << " GeV"
0273                 << ", eion " << syntheticIonization << " GeV"
0274                 << std::endl;
0275     }
0276   }
0277
0278   if (Verbosity() > 0)
0279   {
0280     std::cout << Name() << " " << config.name
0281               << ": inspected " << inspected
0282               << ", modified " << modified
0283               << ", missing particle " << missingParticle
0284               << ", invalid path " << invalidPath
0285               << std::endl;
0286   }
0287 }