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