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   return;
0039 }
0040
0041 double
0042 PHG4SteppingAction::GetScintLightYield(const G4Step* step)
0043 {
0044   //      pirated from G4Scintillation::PostStepDoIt()
0045
0046   double light_yield = 0;
0047
0048   const G4Track* aTrack = step->GetTrack();
0049   const G4Material* aMaterial = aTrack->GetMaterial();
0050   G4MaterialPropertiesTable* aMaterialPropertiesTable =
0051       aMaterial->GetMaterialPropertiesTable();
0052   if (!aMaterialPropertiesTable)
0053   {
0054     const std::string& mname(aMaterial->GetName());
0055
0056     std::set<std::string>::const_iterator it =
0057         m_ScintLightYieldMissingMaterialSet.find(mname);
0058
0059     if (it == m_ScintLightYieldMissingMaterialSet.end())
0060     {
0061       m_ScintLightYieldMissingMaterialSet.insert(mname);
0062
0063       std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
0064                 << "can not find Material Properties Table for material " << mname
0065                 << ", will assume it do NOT scintillate. "
0066                 << "Please ignore this warning if you do not expect scintillation light from "
0067                 << mname << std::endl;
0068     }
0069
0070     return 0.;
0071   }
0072
0073   if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
0074   {
0075     light_yield = aMaterialPropertiesTable->GetConstProperty("SCINTILLATIONYIELD") * GetVisibleEnergyDeposition(step) * GeV;
0076
0077     return light_yield;
0078   }  //  if (aMaterialPropertiesTable->ConstPropertyExists("SCINTILLATIONYIELD"))
0079   const std::string& mname(aMaterial->GetName());
0080
0081   std::set<std::string>::const_iterator it =
0082       m_ScintLightYieldMissingMaterialSet.find(mname);
0083
0084   if (it == m_ScintLightYieldMissingMaterialSet.end())
0085   {
0086     m_ScintLightYieldMissingMaterialSet.insert(mname);
0087
0088     std::cout << "PHG4SteppingAction::GetScintLightYield - WARNING - "
0089               << "can not find scintillation light yield for material " << mname
0090               << ", will assume it do NOT scintillate. "
0091               << "Please ignore this warning if you do not expect scintillation light from "
0092               << mname << std::endl;
0093   }
0094
0095   return light_yield;  // whatever was set as default at beginning
0096 }
0097
0098 double
0099 PHG4SteppingAction::GetVisibleEnergyDeposition(const G4Step* step)
0100 {
0101   G4EmSaturation* emSaturation = G4LossTableManager::Instance()->EmSaturation();
0102   if (emSaturation)
0103   {
0104     if (m_Verbosity)
0105     {
0106       emSaturation->SetVerbose(m_Verbosity);
0107     }
0108     double visen = emSaturation->VisibleEnergyDepositionAtAStep(step) / GeV;
0109     return visen;
0110   }
0111
0112   std::cout
0113       << "PHG4SteppingAction::GetScintLightYield - ERROR - can NOT initialize G4EmSaturation!"
0114       << std::endl;
0115
0116   return 0.;
0117 }
0118
0119 void PHG4SteppingAction::StoreLocalCoordinate(PHG4Hit* hit, const G4Step* aStep,
0120                                               const bool do_prepoint, const bool do_postpoint)
0121 {
0122   assert(hit);
0123   assert(aStep);
0124
0125   G4StepPoint* preStepPoint = aStep->GetPreStepPoint();
0126   const G4TouchableHandle& theTouchable = preStepPoint->GetTouchableHandle();
0127
0128   if (do_prepoint)
0129   {
0130     const G4ThreeVector& worldPosition = preStepPoint->GetPosition();
0131     G4ThreeVector localPosition =
0132         theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0133
0134     hit->set_local_x(0, localPosition.x() / cm);
0135     hit->set_local_y(0, localPosition.y() / cm);
0136     hit->set_local_z(0, localPosition.z() / cm);
0137   }
0138   if (do_postpoint)
0139   {
0140     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0141
0142     const G4ThreeVector& worldPosition = postPoint->GetPosition();
0143     G4ThreeVector localPosition =
0144         theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
0145
0146     hit->set_local_x(1, localPosition.x() / cm);
0147     hit->set_local_y(1, localPosition.y() / cm);
0148     hit->set_local_z(1, localPosition.z() / cm);
0149   }
0150 }
0151
0152 void PHG4SteppingAction::SetLightCorrection(const double inner_radius, const double inner_corr, const double outer_radius, const double outer_corr)
0153 {
0154   m_LightBalanceInnerRadius = inner_radius;
0155   m_LightBalanceInnerCorr = inner_corr;
0156   m_LightBalanceOuterRadius = outer_radius;
0157   m_LightBalanceOuterCorr = outer_corr;
0158   return;
0159 }
0160
0161 double PHG4SteppingAction::GetLightCorrection(const double xpos, const double ypos) const
0162 {
0163   double r = sqrt(xpos * xpos + ypos * ypos);
0164   double correction = GetLightCorrection(r);
0165   return correction;
0166 }
0167
0168 double PHG4SteppingAction::GetLightCorrection(const double r) const
0169 {
0170   double correction = 1.;
0171   if (ValidCorrection())
0172   {
0173     double slope = (m_LightBalanceOuterCorr - m_LightBalanceInnerCorr) / (m_LightBalanceOuterRadius - m_LightBalanceInnerRadius);
0174     double b = m_LightBalanceInnerCorr - m * m_LightBalanceInnerRadius;
0175     correction = slope * r + b;
0176     correction = std::min(correction, 1.);
0177     correction = std::max(correction, 0.);
0178   }
0179   return correction;
0180 }
0181
0182 bool PHG4SteppingAction::ValidCorrection() const
0183 {
0184   if (std::isfinite(m_LightBalanceOuterRadius) &&
0185       std::isfinite(m_LightBalanceInnerRadius) &&
0186       std::isfinite(m_LightBalanceOuterCorr) &&
0187       std::isfinite(m_LightBalanceInnerCorr))
0188   {
0189     return true;
0190   }
0191   return false;
0192 }