g4simulation/g4main/PHG4SteppingAction.cc

0001 /*!
0002  * \file PHG4SteppingAction.cc
0003  * \brief
0004  * \author Jin Huang <jhuang@bnl.gov>
0005  * \version $Revision: 1.1 $
0006  * \date $Date: 2015/01/07 23:50:05 $
0007  */
0008
0009 #include "PHG4SteppingAction.h"
0010
0011 #include "PHG4Hit.h"
0012
0013 #include <Geant4/G4AffineTransform.hh>  // for G4AffineTransform
0014 #include <Geant4/G4EmSaturation.hh>
0015 #include <Geant4/G4LossTableManager.hh>
0016 #include <Geant4/G4Material.hh>
0017 #include <Geant4/G4MaterialPropertiesTable.hh>  // for G4MaterialProperties...
0018 #include <Geant4/G4NavigationHistory.hh>
0019 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHa...
0020 #include <Geant4/G4Step.hh>
0021 #include <Geant4/G4StepPoint.hh>
0022 #include <Geant4/G4String.hh>  // for G4String
0023 #include <Geant4/G4SystemOfUnits.hh>
0024 #include <Geant4/G4ThreeVector.hh>
0025 #include <Geant4/G4TouchableHandle.hh>  // for G4TouchableHandle
0026 #include <Geant4/G4Track.hh>
0027 #include <Geant4/G4VTouchable.hh>  // for G4VTouchable
0028
0029 #include <algorithm>
0030 #include <cassert>
0031 #include <cmath>  // for isfinite, NAN, sqrt
0032 #include <iostream>
0033
0034 PHG4SteppingAction::PHG4SteppingAction(const std::string& name, const int i)
0035   : m_Verbosity(i)
0036   , m_Name(name)
0037 {
0038 }
0039
0040 double
0041 PHG4SteppingAction::GetScintLightYield(const G4Step* step)
0042 {
0043   //      pirated from G4Scintillation::PostStepDoIt()
0044
0045   double light_yield = 0;
0046
0047   const G4Track* aTrack = step->GetTrack();
0048   const G4Material* aMaterial = aTrack->GetMaterial();
0049   G4MaterialPropertiesTable* aMaterialPropertiesTable =
0050       aMaterial->GetMaterialPropertiesTable();
0051   if (!aMaterialPropertiesTable)
0052   {
0053     const std::string& mname(aMaterial->GetName());
0054
0055     std::set<std::string>::const_iterator it =
0056         m_ScintLightYieldMissingMaterialSet.find(mname);
0057
0058     if (it == m_ScintLightYieldMissingMaterialSet.end())
0059     {
0060       m_ScintLightYieldMissingMaterialSet.insert(mname);
0061
0062       std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
0063                 << "can not find Material Properties Table for material " << mname
0064                 << ", will assume it do NOT scintillate. "
0065                 << "Please ignore this warning if you do not expect scintillation light from "
0066                 << mname << std::endl;
0067     }
0068
0069     return 0.;
0070   }
0071
0072   if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
0073   {
0074     light_yield = aMaterialPropertiesTable->GetConstProperty("SCINTILLATIONYIELD") * GetVisibleEnergyDeposition(step) * GeV;
0075
0076     return light_yield;
0077   }  //  if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
0078   else
0079   {
0080     const std::string& mname(aMaterial->GetName());
0081
0082     std::set<std::string>::const_iterator it =
0083         m_ScintLightYieldMissingMaterialSet.find(mname);
0084
0085     if (it == m_ScintLightYieldMissingMaterialSet.end())
0086     {
0087       m_ScintLightYieldMissingMaterialSet.insert(mname);
0088
0089       std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
0090                 << "can not find scintillation light yield for material " << mname
0091                 << ", will assume it do NOT scintillate. "
0092                 << "Please ignore this warning if you do not expect scintillation light from "
0093                 << mname << std::endl;
0094     }
0095
0096     return 0.;
0097   }
0098
0099   return light_yield;
0100 }
0101
0102 double
0103 PHG4SteppingAction::GetVisibleEnergyDeposition(const G4Step* step)
0104 {
0105   G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
0106   if (emSaturation)
0107   {
0108     if (m_Verbosity)
0109     {
0110       emSaturation->SetVerbose(m_Verbosity);
0111     }
0112     double visen = emSaturation->VisibleEnergyDepositionAtAStep(step) / GeV;
0113     return visen;
0114   }
0115   else
0116   {
0117     std::cout
0118         << "PHG4SteppingAction::GetScintLightYield - ERROR - can NOT initialize G4EmSaturation!"
0119         << std::endl;
0120
0121     return 0.;
0122   }
0123 }
0124
0125 void PHG4SteppingAction::StoreLocalCoordinate(PHG4Hit* hit, const G4Step* aStep,
0126                                               const bool do_prepoint, const bool do_postpoint)
0127 {
0128   assert(hit);
0129   assert(aStep);
0130
0131   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0132   const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
0133
0134   if (do_prepoint)
0135   {
0136     const G4ThreeVector& worldPosition = preStepPoint->GetPosition();
0137     G4ThreeVector localPosition =
0138         theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0139
0140     hit->set_local_x(0, localPosition.x() / cm);
0141     hit->set_local_y(0, localPosition.y() / cm);
0142     hit->set_local_z(0, localPosition.z() / cm);
0143   }
0144   if (do_postpoint)
0145   {
0146     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0147
0148     const G4ThreeVector& worldPosition = postPoint->GetPosition();
0149     G4ThreeVector localPosition =
0150         theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0151
0152     hit->set_local_x(1, localPosition.x() / cm);
0153     hit->set_local_y(1, localPosition.y() / cm);
0154     hit->set_local_z(1, localPosition.z() / cm);
0155   }
0156 }
0157
0158 void PHG4SteppingAction::SetLightCorrection(const double inner_radius, const double inner_corr, const double outer_radius, const double outer_corr)
0159 {
0160   m_LightBalanceInnerRadius = inner_radius;
0161   m_LightBalanceInnerCorr = inner_corr;
0162   m_LightBalanceOuterRadius = outer_radius;
0163   m_LightBalanceOuterCorr = outer_corr;
0164   return;
0165 }
0166
0167 double PHG4SteppingAction::GetLightCorrection(const double xpos, const double ypos) const
0168 {
0169   double r = sqrt(xpos * xpos + ypos * ypos);
0170   double correction = GetLightCorrection(r);
0171   return correction;
0172 }
0173
0174 double PHG4SteppingAction::GetLightCorrection(const double r) const
0175 {
0176   double correction = 1.;
0177   if (ValidCorrection())
0178   {
0179     double slope = (m_LightBalanceOuterCorr - m_LightBalanceInnerCorr) / (m_LightBalanceOuterRadius - m_LightBalanceInnerRadius);
0180     double b = m_LightBalanceInnerCorr - m * m_LightBalanceInnerRadius;
0181     correction = slope * r + b;
0182     correction = std::min(correction, 1.);
0183     correction = std::max(correction, 0.);
0184   }
0185   return correction;
0186 }
0187
0188 bool PHG4SteppingAction::ValidCorrection() const
0189 {
0190   if (std::isfinite(m_LightBalanceOuterRadius) &&
0191       std::isfinite(m_LightBalanceInnerRadius) &&
0192       std::isfinite(m_LightBalanceOuterCorr) &&
0193       std::isfinite(m_LightBalanceInnerCorr))
0194   {
0195     return true;
0196   }
0197   return false;
0198 }