g4simulation/g4epd/PHG4EPDSteppingAction.cc

0001 /* vim: set sw=2 ft=cpp: */
0002
0003 #include "PHG4EPDSteppingAction.h"
0004
0005 #include "PHG4EPDDetector.h"
0006
0007 #include <g4detectors/PHG4StepStatusDecode.h>
0008
0009 #include <phool/getClass.h>
0010
0011 #include <g4main/PHG4Hit.h>
0012 #include <g4main/PHG4HitContainer.h>
0013 #include <g4main/PHG4Hitv1.h>
0014 #include <g4main/PHG4Shower.h>
0015 #include <g4main/PHG4SteppingAction.h>
0016 #include <g4main/PHG4TrackUserInfoV1.h>
0017
0018 #include <TSystem.h>
0019
0020 #include <Geant4/G4ParticleDefinition.hh>
0021 #include <Geant4/G4Step.hh>
0022 #include <Geant4/G4StepPoint.hh>
0023 #include <Geant4/G4StepStatus.hh>
0024 #include <Geant4/G4String.hh>  // for G4String
0025 #include <Geant4/G4SystemOfUnits.hh>
0026 #include <Geant4/G4ThreeVector.hh>
0027 #include <Geant4/G4TouchableHandle.hh>
0028 #include <Geant4/G4Track.hh>
0029 #include <Geant4/G4TrackStatus.hh>
0030 #include <Geant4/G4Types.hh>
0031 #include <Geant4/G4VTouchable.hh>
0032 #include <Geant4/G4VUserTrackInformation.hh>
0033
0034 #include <cstdint>  // for int32_t
0035 #include <iostream>
0036 #include <string>
0037
0038 class G4VPhysicalVolume;
0039
0040 PHG4EPDSteppingAction::PHG4EPDSteppingAction(PHG4EPDDetector* detector,
0041                                              const PHParameters* /*unused*/)
0042   : PHG4SteppingAction(detector->GetName())
0043   , m_Detector(detector)
0044 {
0045 }
0046
0047 PHG4EPDSteppingAction::~PHG4EPDSteppingAction()
0048 {
0049   delete m_Hit;
0050 }
0051
0052 bool PHG4EPDSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0053 {
0054   G4StepPoint* prestep = aStep->GetPreStepPoint();
0055   G4TouchableHandle prehandle = prestep->GetTouchableHandle();
0056
0057   G4VPhysicalVolume* volume = prehandle->GetVolume();
0058
0059   // m_Detector->IsInDetector(volume)
0060   // returns
0061   //  0 is outside of EPD
0062   //  1 is inside scintillator
0063   // -1 is inside support structure (dead material)
0064   int whichactive = m_Detector->IsInDetector(volume);
0065   if (!whichactive)
0066   {
0067     return false;
0068   }
0069
0070   G4double deposit = aStep->GetTotalEnergyDeposit() / GeV;
0071   G4double ionising = deposit - aStep->GetNonIonizingEnergyDeposit() / GeV;
0072   G4double light_yield = GetVisibleEnergyDeposition(aStep) / GeV;
0073
0074   //  G4StepStatus prestatus = prestep->GetStepStatus();
0075
0076   int32_t tile_id = m_Detector->module_id_for(volume);
0077
0078   G4Track const* track = aStep->GetTrack();
0079
0080   G4ParticleDefinition const* particle = track->GetParticleDefinition();
0081
0082   bool geantino = (particle->GetPDGEncoding() == 0 && particle->GetParticleName().find("geantino") != std::string::npos);
0083
0084   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0085   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0086
0087   switch (prePoint->GetStepStatus())
0088   {
0089   case fPostStepDoItProc:
0090     if (m_SavePostStepStatus != fGeomBoundary)
0091     {
0092       break;
0093     }
0094     else
0095     {
0096       std::cout << GetName() << ": New Hit for  " << std::endl;
0097       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0098                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0099                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0100     }
0101     [[fallthrough]];
0102   case fGeomBoundary:
0103   case fUndefined:
0104     if (m_Hit == nullptr)
0105     {
0106       m_Hit = new PHG4Hitv1();
0107     }
0108
0109     // only for active columes (scintillators)
0110     if (whichactive > 0)
0111     {
0112       m_Hit->set_eion(0);
0113       m_Hit->set_light_yield(0);
0114     }
0115     m_Hit->set_x(0, prestep->GetPosition().x() / cm);
0116     m_Hit->set_y(0, prestep->GetPosition().y() / cm);
0117     m_Hit->set_z(0, prestep->GetPosition().z() / cm);
0118     m_Hit->set_t(0, prestep->GetGlobalTime() / nanosecond);
0119
0120     m_Hit->set_trkid(track->GetTrackID());
0121
0122     if (PHG4TrackUserInfoV1* userinfo = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation()))
0123     {
0124       m_Hit->set_trkid(userinfo->GetUserTrackId());
0125
0126       userinfo->GetShower()->add_g4hit_id(m_HitContainer->GetID(), m_Hit->get_hit_id());
0127     }
0128
0129     m_Hit->set_edep(0);
0130     break;
0131   default:
0132     break;
0133   }
0134
0135   G4StepPoint* poststep = aStep->GetPostStepPoint();
0136   const G4ThreeVector postpos = poststep->GetPosition();
0137   m_SavePostStepStatus = postPoint->GetStepStatus();
0138   m_Hit->set_edep(m_Hit->get_edep() + deposit);
0139   if (whichactive > 0)
0140   {
0141     m_Hit->set_eion(m_Hit->get_eion() + ionising);
0142     m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield * GetLightCorrection(postpos.x(), postpos.y()));
0143   }
0144
0145   if (postPoint->GetStepStatus() != fGeomBoundary &&
0146       postPoint->GetStepStatus() != fWorldBoundary &&
0147       postPoint->GetStepStatus() != fAtRestDoItProc &&
0148       track->GetTrackStatus() != fStopAndKill)
0149   {
0150     return true;
0151   }
0152   if (m_Hit->get_edep() <= 0 && !geantino)
0153   {
0154     m_Hit->Reset();
0155
0156     return true;
0157   }
0158
0159   m_Hit->set_x(1, poststep->GetPosition().x() / cm);
0160   m_Hit->set_y(1, poststep->GetPosition().y() / cm);
0161   m_Hit->set_z(1, poststep->GetPosition().z() / cm);
0162   m_Hit->set_t(1, poststep->GetGlobalTime() / nanosecond);
0163
0164   PHG4TrackUserInfoV1* userinfo = dynamic_cast<PHG4TrackUserInfoV1*>(track->GetUserInformation());
0165
0166   if (userinfo != nullptr)
0167   {
0168     userinfo->SetKeep(1);
0169   }
0170
0171   if (geantino)
0172   {
0173     m_Hit->set_edep(-1.);
0174     m_Hit->set_eion(-1.);
0175     m_Hit->set_light_yield(-1.);
0176   }
0177
0178   m_HitContainer->AddHit(tile_id, m_Hit);
0179
0180   m_Hit = nullptr;
0181
0182   return true;
0183 }
0184
0185 void PHG4EPDSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0186 {
0187   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0188
0189   m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
0190   // if we do not find the node it's messed up.
0191   if (!m_HitContainer)
0192   {
0193     std::cout << "PHG4EPDSteppingAction::SetTopNode - unable to find hit node " << m_HitNodeName << std::endl;
0194     gSystem->Exit(1);
0195   }
0196   // this is perfectly fine if support hits are disabled
0197   if (!m_SupportHitContainer)
0198   {
0199     if (Verbosity() > 0)
0200     {
0201       std::cout << "PHG4EPDSteppingAction::SetTopNode - unable to find support hit node " << m_SupportNodeName << std::endl;
0202     }
0203   }
0204 }
0205
0206 void PHG4EPDSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
0207 {
0208   if (type == "G4HIT")
0209   {
0210     m_HitNodeName = name;
0211     return;
0212   }
0213   if (type == "G4HIT_SUPPORT")
0214   {
0215     m_SupportNodeName = name;
0216     return;
0217   }
0218   std::cout << "Invalid output hit node type " << type << std::endl;
0219   gSystem->Exit(1);
0220   return;
0221 }