g4simulation/g4detectors/PHG4BlockSteppingAction.cc

0001 #include "PHG4BlockSteppingAction.h"
0002
0003 #include "PHG4BlockDetector.h"
0004 #include "PHG4StepStatusDecode.h"
0005
0006 #include <phparameter/PHParameters.h>
0007
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010 #include <g4main/PHG4Hitv1.h>
0011 #include <g4main/PHG4Shower.h>
0012 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0013 #include <g4main/PHG4TrackUserInfoV1.h>
0014
0015 #include <phool/getClass.h>
0016
0017 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0018 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0019 #include <Geant4/G4Step.hh>
0020 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0021 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fPostSt...
0022 #include <Geant4/G4String.hh>      // for G4String
0023 #include <Geant4/G4SystemOfUnits.hh>
0024 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0025 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0026 #include <Geant4/G4Track.hh>                  // for G4Track
0027 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0028 #include <Geant4/G4Types.hh>                  // for G4double
0029 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0030 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0031 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0032
0033 #include <cmath>    // for isfinite
0034 #include <cstdlib>  // for exit
0035 #include <iostream>
0036 #include <string>  // for operator<<, string
0037
0038 class PHCompositeNode;
0039
0040 using namespace std;
0041 //____________________________________________________________________________..
0042 PHG4BlockSteppingAction::PHG4BlockSteppingAction(PHG4BlockDetector* detector, const PHParameters* parameters)
0043   : PHG4SteppingAction(detector->GetName())
0044   , m_Detector(detector)
0045   , m_Params(parameters)
0046   , m_HitContainer(nullptr)
0047   , m_Hit(nullptr)
0048   , m_SaveShower(nullptr)
0049   , m_SaveVolPre(nullptr)
0050   , m_SaveVolPost(nullptr)
0051   , m_SaveTrackId(-1)
0052   , m_SavePreStepStatus(-1)
0053   , m_SavePostStepStatus(-1)
0054   , m_ActiveFlag(m_Params->get_int_param("active"))
0055   , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
0056   , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
0057 {
0058 }
0059
0060 PHG4BlockSteppingAction::~PHG4BlockSteppingAction()
0061 {
0062   // if the last hit was a zero energie deposit hit, it is just reset
0063   // and the memory is still allocated, so we need to delete it here
0064   // if the last hit was saved, hit is a nullptr pointer which are
0065   // legal to delete (it results in a no operation)
0066   delete m_Hit;
0067 }
0068
0069 //____________________________________________________________________________..
0070 bool PHG4BlockSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0071 {
0072   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0073   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0074   // get volume of the current step
0075   G4VPhysicalVolume* volume = touch->GetVolume();
0076
0077   if (!m_Detector->IsInBlock(volume))
0078   {
0079     return false;
0080   }
0081
0082   // collect energy and track length step by step
0083   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0084   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0085   const G4Track* aTrack = aStep->GetTrack();
0086
0087   // if this block stops everything, just put all kinetic energy into edep
0088   if (m_BlackHoleFlag)
0089   {
0090     edep = aTrack->GetKineticEnergy() / GeV;
0091     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0092     killtrack->SetTrackStatus(fStopAndKill);
0093   }
0094
0095   // make sure we are in a volume
0096   if (m_ActiveFlag)
0097   {
0098     int layer_id = m_Detector->get_Layer();
0099     bool geantino = false;
0100     // the check for the pdg code speeds things up, I do not want to make
0101     // an expensive string compare for every track when we know
0102     // geantino or chargedgeantino has pid=0
0103     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0104         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
0105     {
0106       geantino = true;
0107     }
0108     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0109     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0110     //       cout << "track id " << aTrack->GetTrackID() << endl;
0111     //       cout << "time prepoint: " << prePoint->GetGlobalTime() << endl;
0112     //       cout << "time postpoint: " << postPoint->GetGlobalTime() << endl;
0113     int prepointstatus = prePoint->GetStepStatus();
0114     if (prepointstatus == fGeomBoundary ||
0115         prepointstatus == fUndefined ||
0116         (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
0117         m_UseG4StepsFlag > 0)
0118     {
0119       if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
0120       {
0121         cout << GetName() << ": New Hit for  " << endl;
0122         cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0123              << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0124              << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0125              << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
0126         cout << "last track: " << m_SaveTrackId
0127              << ", current trackid: " << aTrack->GetTrackID() << endl;
0128         cout << "phys pre vol: " << volume->GetName()
0129              << " post vol : " << touchpost->GetVolume()->GetName() << endl;
0130         cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0131              << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
0132       }
0133       if (!m_Hit)
0134       {
0135         m_Hit = new PHG4Hitv1();
0136       }
0137       m_Hit->set_layer(layer_id);
0138       // here we set the entrance values in cm
0139       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0140       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0141       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0142       // time in ns
0143       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0144       // set the track ID
0145       m_Hit->set_trkid(aTrack->GetTrackID());
0146       m_SaveTrackId = aTrack->GetTrackID();
0147
0148       // set the initial energy deposit
0149       m_Hit->set_edep(0);
0150       if (!geantino && !m_BlackHoleFlag && m_ActiveFlag)
0151       {
0152         m_Hit->set_eion(0);
0153       }
0154       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0155       {
0156         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0157         {
0158           m_Hit->set_trkid(pp->GetUserTrackId());
0159           m_Hit->set_shower_id(pp->GetShower()->get_id());
0160           m_SaveShower = pp->GetShower();
0161         }
0162       }
0163     }
0164
0165     // some sanity checks for inconsistencies
0166     // check if this hit was created, if not print out last post step status
0167     if (!m_Hit || !isfinite(m_Hit->get_x(0)))
0168     {
0169       cout << GetName() << ": hit was not created" << endl;
0170       cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0171            << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0172            << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0173            << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << endl;
0174       cout << "last track: " << m_SaveTrackId
0175            << ", current trackid: " << aTrack->GetTrackID() << endl;
0176       cout << "phys pre vol: " << volume->GetName()
0177            << " post vol : " << touchpost->GetVolume()->GetName() << endl;
0178       cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0179            << " previous phys post vol: " << m_SaveVolPost->GetName() << endl;
0180       exit(1);
0181     }
0182     m_SavePostStepStatus = postPoint->GetStepStatus();
0183     // check if track id matches the initial one when the hit was created
0184     if (aTrack->GetTrackID() != m_SaveTrackId)
0185     {
0186       cout << "hits do not belong to the same track" << endl;
0187       cout << "saved track: " << m_SaveTrackId
0188            << ", current trackid: " << aTrack->GetTrackID()
0189            << endl;
0190       exit(1);
0191     }
0192     m_SavePreStepStatus = prePoint->GetStepStatus();
0193     m_SavePostStepStatus = postPoint->GetStepStatus();
0194     m_SaveVolPre = volume;
0195     m_SaveVolPost = touchpost->GetVolume();
0196
0197     // here we just update the exit values, it will be overwritten
0198     // for every step until we leave the volume or the particle
0199     // ceases to exist
0200     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0201     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0202     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0203
0204     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0205     // sum up the energy to get total deposited
0206     m_Hit->set_edep(m_Hit->get_edep() + edep);
0207     // update ionization energy only for active volumes, not for black holes or geantinos
0208     // if the hit is created without eion, get_eion() returns NAN
0209     // if you see NANs check the creation of the hit
0210     if (m_ActiveFlag && !m_BlackHoleFlag && !geantino)
0211     {
0212       m_Hit->set_eion(m_Hit->get_eion() + eion);
0213     }
0214     if (geantino)
0215     {
0216       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0217     }
0218     if (edep > 0)
0219     {
0220       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0221       {
0222         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0223         {
0224           pp->SetKeep(1);  // we want to keep the track
0225         }
0226       }
0227     }
0228     // if any of these conditions is true this is the last step in
0229     // this volume and we need to save the hit
0230     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0231     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0232     // (happens when your detector goes outside world volume)
0233     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0234     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0235     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0236     if (postPoint->GetStepStatus() == fGeomBoundary ||
0237         postPoint->GetStepStatus() == fWorldBoundary ||
0238         postPoint->GetStepStatus() == fAtRestDoItProc ||
0239         aTrack->GetTrackStatus() == fStopAndKill ||
0240         m_UseG4StepsFlag > 0)
0241     {
0242       // save only hits with energy deposit (or -1 for geantino)
0243       if (m_Hit->get_edep())
0244       {
0245         m_HitContainer->AddHit(layer_id, m_Hit);
0246         if (m_SaveShower)
0247         {
0248           m_SaveShower->add_g4hit_id(m_HitContainer->GetID(), m_Hit->get_hit_id());
0249         }
0250         // ownership has been transferred to container, set to null
0251         // so we will create a new hit for the next track
0252         m_Hit = nullptr;
0253       }
0254       else
0255       {
0256         // if this hit has no energy deposit, just reset it for reuse
0257         // this means we have to delete it in the dtor. If this was
0258         // the last hit we processed the memory is still allocated
0259         m_Hit->Reset();
0260       }
0261     }
0262     // return true to indicate the hit was used
0263     return true;
0264   }
0265   else
0266   {
0267     return false;
0268   }
0269 }
0270
0271 //____________________________________________________________________________..
0272 void PHG4BlockSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0273 {
0274   string hitnodename;
0275   if (m_Detector->SuperDetector() != "NONE")
0276   {
0277     hitnodename = "G4HIT_" + m_Detector->SuperDetector();
0278   }
0279   else
0280   {
0281     hitnodename = "G4HIT_" + m_Detector->GetName();
0282   }
0283
0284   // now look for the map and grab a pointer to it.
0285   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0286
0287   // if we do not find the node we need to make it.
0288   if (!m_HitContainer)
0289   {
0290     std::cout << "PHG4BlockSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0291   }
0292 }