g4simulation/g4tpc/PHG4TpcEndCapSteppingAction.cc

0001 #include "PHG4TpcEndCapSteppingAction.h"
0002
0003 #include "PHG4TpcEndCapDetector.h"
0004
0005 #include <phparameter/PHParameters.h>
0006
0007 #include <g4detectors/PHG4StepStatusDecode.h>
0008
0009 #include <g4main/PHG4Hit.h>
0010 #include <g4main/PHG4HitContainer.h>
0011 #include <g4main/PHG4Hitv1.h>
0012 #include <g4main/PHG4Shower.h>
0013 #include <g4main/PHG4SteppingAction.h>
0014 #include <g4main/PHG4TrackUserInfoV1.h>
0015
0016 #include <phool/getClass.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>
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/G4VPhysicalVolume.hh>
0033 #include <Geant4/G4VTouchable.hh>
0034 #include <Geant4/G4VUserTrackInformation.hh>
0035
0036 #include <cmath>
0037 #include <iostream>
0038 #include <string>
0039
0040 class PHCompositeNode;
0041
0042 //____________________________________________________________________________..
0043 PHG4TpcEndCapSteppingAction::PHG4TpcEndCapSteppingAction(PHG4TpcEndCapDetector *detector, const PHParameters *parameters)
0044   : PHG4SteppingAction(detector->GetName())
0045   , m_Detector(detector)
0046   , m_Params(parameters)
0047   , m_ActiveFlag(m_Params->get_int_param("active"))
0048   , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
0049 {
0050 }
0051
0052 //____________________________________________________________________________..
0053 PHG4TpcEndCapSteppingAction::~PHG4TpcEndCapSteppingAction()
0054 {
0055   // if the last hit was a zero energie deposit hit, it is just reset
0056   // and the memory is still allocated, so we need to delete it here
0057   // if the last hit was saved, hit is a nullptr pointer which are
0058   // legal to delete (it results in a no operation)
0059   delete m_Hit;
0060 }
0061
0062 //____________________________________________________________________________..
0063 // This is the implementation of the G4 UserSteppingAction
0064 bool PHG4TpcEndCapSteppingAction::UserSteppingAction(const G4Step *aStep, bool /*was_used*/)
0065 {
0066   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0067   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0068   // get volume of the current step
0069   G4VPhysicalVolume *volume = touch->GetVolume();
0070   // IsInDetector(volume) returns
0071   //  == 0 outside of detector
0072   //   > 0 for hits in active volume
0073   //  < 0 for hits in passive material
0074   int whichactive = m_Detector->IsInDetector(volume);
0075   if (!whichactive)
0076   {
0077     return false;
0078   }
0079   // collect energy and track length step by step
0080   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0081   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0082   const G4Track *aTrack = aStep->GetTrack();
0083   // if this detector stops everything, just put all kinetic energy into edep
0084   if (m_BlackHoleFlag)
0085   {
0086     edep = aTrack->GetKineticEnergy() / GeV;
0087     G4Track *killtrack = const_cast<G4Track *>(aTrack);
0088     killtrack->SetTrackStatus(fStopAndKill);
0089   }
0090   // we use here only one detector in this simple example
0091   // if you deal with multiple detectors in this stepping action
0092   // the detector id can be used to distinguish between them
0093   // hits can easily be analyzed later according to their detector id
0094   int detector_id = 0;  // we use here only one detector in this simple example
0095   bool geantino = false;
0096   // the check for the pdg code speeds things up, I do not want to make
0097   // an expensive string compare for every track when we know
0098   // geantino or chargedgeantino has pid=0
0099   if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0100       aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
0101           std::string::npos)  // this also accounts for "chargedgeantino"
0102   {
0103     geantino = true;
0104   }
0105   G4StepPoint *prePoint = aStep->GetPreStepPoint();
0106   G4StepPoint *postPoint = aStep->GetPostStepPoint();
0107
0108   // Here we have to decide if we need to create a new hit.  Normally this should
0109   // only be neccessary if a G4 Track enters a new volume or is freshly created
0110   // For this we look at the step status of the prePoint (beginning of the G4 Step).
0111   // This should be either fGeomBoundary (G4 Track crosses into volume) or
0112   // fUndefined (G4 Track newly created)
0113   // Sadly over the years with different G4 versions we have observed cases where
0114   // G4 produces "impossible hits" which we try to catch here
0115   // These errors were always rare and it is not clear if they still exist but we
0116   // still check for them for safety. We can reproduce G4 runs identically (if given
0117   // the sequence of random number seeds you find in the log), the printouts help
0118   // us giving the G4 support information about those failures
0119   //
0120   switch (prePoint->GetStepStatus())
0121   {
0122   case fPostStepDoItProc:
0123     if (m_SavePostStepStatus != fGeomBoundary)
0124     {
0125       // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
0126       // a new volume, just proceed here
0127       break;
0128     }
0129     else
0130     {
0131       // this is an impossible G4 Step print out diagnostic to help debug, not sure if
0132       // this is still with us
0133       std::cout << GetName() << ": New Hit for  " << std::endl;
0134       std::cout << "prestep status: "
0135                 << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0136                 << ", poststep status: "
0137                 << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0138                 << ", last pre step status: "
0139                 << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0140                 << ", last post step status: "
0141                 << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0142       std::cout << "last track: " << m_SaveTrackId
0143                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0144       std::cout << "phys pre vol: " << volume->GetName()
0145                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0146       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0147                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0148     }
0149     [[fallthrough]];
0150     // These are the normal cases
0151   case fGeomBoundary:
0152   case fUndefined:
0153     if (!m_Hit)
0154     {
0155       m_Hit = new PHG4Hitv1();
0156     }
0157     m_Hit->set_layer(detector_id);
0158     // here we set the entrance values in cm
0159     m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0160     m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0161     m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0162     // time in ns
0163     m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0164     // set the track ID
0165     m_Hit->set_trkid(aTrack->GetTrackID());
0166     m_SaveTrackId = aTrack->GetTrackID();
0167     // set the initial energy deposit
0168     m_EdepSum = 0;
0169     // implement your own here://
0170     // add the properties you are interested in via set_XXX methods
0171     // you can find existing set methods in $OFFLINE_MAIN/include/g4main/PHG4Hit.h
0172     // this is initialization of your value. This is not needed you can just set the final
0173     // value at the last step in this volume later one
0174     if (whichactive > 0)
0175     {
0176       m_EionSum = 0;  // assuming the ionization energy is only needed for active
0177                       // volumes (scintillators)
0178       m_Hit->set_eion(0);
0179       m_SaveHitContainer = m_HitContainer;
0180     }
0181     else
0182     {
0183       std::cout << "implement stuff for whichactive < 0 (inactive volumes)" << std::endl;
0184       gSystem->Exit(1);
0185     }
0186     // this is for the tracking of the truth info
0187     if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0188     {
0189       if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0190       {
0191         m_Hit->set_trkid(pp->GetUserTrackId());
0192         pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0193       }
0194     }
0195     break;
0196   default:
0197     break;
0198   }
0199
0200   // This section is called for every step
0201   // some sanity checks for inconsistencies (aka bugs) we have seen over the years
0202   // check if this hit was created, if not print out last post step status
0203   if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0204   {
0205     std::cout << GetName() << ": hit was not created" << std::endl;
0206     std::cout << "prestep status: "
0207               << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0208               << ", poststep status: "
0209               << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0210               << ", last pre step status: "
0211               << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0212               << ", last post step status: "
0213               << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0214     std::cout << "last track: " << m_SaveTrackId
0215               << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0216     std::cout << "phys pre vol: " << volume->GetName()
0217               << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0218     std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0219               << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0220     // This is fatal - a hit from nowhere. This needs to be looked at and fixed
0221     gSystem->Exit(1);
0222   }
0223   // check if track id matches the initial one when the hit was created
0224   if (aTrack->GetTrackID() != m_SaveTrackId)
0225   {
0226     std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0227     std::cout << "saved track: " << m_SaveTrackId
0228               << ", current trackid: " << aTrack->GetTrackID()
0229               << ", prestep status: " << prePoint->GetStepStatus()
0230               << ", previous post step status: " << m_SavePostStepStatus << std::endl;
0231     // This is fatal - a hit from nowhere. This needs to be looked at and fixed
0232     gSystem->Exit(1);
0233   }
0234
0235   // We need to cache a few things from one step to the next
0236   // to identify impossible hits and subsequent debugging printout
0237   m_SavePreStepStatus = prePoint->GetStepStatus();
0238   m_SavePostStepStatus = postPoint->GetStepStatus();
0239   m_SaveVolPre = volume;
0240   m_SaveVolPost = touchpost->GetVolume();
0241   // here we just update the exit values, it will be overwritten
0242   // for every step until we leave the volume or the particle
0243   // ceases to exist
0244   // sum up the energy to get total deposited
0245   m_EdepSum += edep;
0246   if (whichactive > 0)
0247   {
0248     m_EionSum += eion;
0249   }
0250   // if any of these conditions is true this is the last step in
0251   // this volume and we need to save the hit
0252   // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0253   // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0254   // (happens when your detector goes outside world volume)
0255   // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0256   // aTrack->GetTrackStatus() == fStopAndKill is also set)
0257   // aTrack->GetTrackStatus() == fStopAndKill: track ends
0258   if (postPoint->GetStepStatus() == fGeomBoundary ||
0259       postPoint->GetStepStatus() == fWorldBoundary ||
0260       postPoint->GetStepStatus() == fAtRestDoItProc ||
0261       aTrack->GetTrackStatus() == fStopAndKill)
0262   {
0263     // save only hits with energy deposit (or geantino)
0264     if (m_EdepSum > 0 || geantino)
0265     {
0266       // update values at exit coordinates and set keep flag
0267       // of track to keep
0268       m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0269       m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0270       m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0271       m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0272       if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0273       {
0274         if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0275         {
0276           pp->SetKeep(1);  // we want to keep the track
0277         }
0278       }
0279       if (geantino)
0280       {
0281         //implement your own here://
0282         // if you want to do something special for geantinos (normally you do not)
0283         m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way
0284                               // geantinos survive the g4hit compression
0285         if (whichactive > 0)
0286         {
0287           m_Hit->set_eion(-1);
0288         }
0289       }
0290       else
0291       {
0292         m_Hit->set_edep(m_EdepSum);
0293       }
0294       //implement your own here://
0295       // what you set here will be saved in the output
0296       if (whichactive > 0)
0297       {
0298         m_Hit->set_eion(m_EionSum);
0299       }
0300       m_SaveHitContainer->AddHit(detector_id, m_Hit);
0301       // ownership has been transferred to container, set to null
0302       // so we will create a new hit for the next track
0303       m_Hit = nullptr;
0304     }
0305     else
0306     {
0307       // if this hit has no energy deposit, just reset it for reuse
0308       // this means we have to delete it in the dtor. If this was
0309       // the last hit we processed the memory is still allocated
0310       m_Hit->Reset();
0311     }
0312   }
0313   // return true to indicate the hit was used
0314   return true;
0315 }
0316
0317 //____________________________________________________________________________..
0318 void PHG4TpcEndCapSteppingAction::SetInterfacePointers(PHCompositeNode *topNode)
0319 {
0320   // now look for the map and grab a pointer to it.
0321   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0322   // if we do not find the node we need to make it.
0323   if (!m_HitContainer)
0324   {
0325     std::cout << "PHG4TpcEndCapSteppingAction::SetTopNode - unable to find "
0326               << m_HitNodeName << std::endl;
0327     gSystem->Exit(1);
0328   }
0329 }
0330
0331 void PHG4TpcEndCapSteppingAction::SetHitNodeName(const std::string &type, const std::string &name)
0332 {
0333   if (type == "G4HIT")
0334   {
0335     m_HitNodeName = name;
0336     return;
0337   }
0338   std::cout << "Invalid output hit node type " << type << std::endl;
0339   gSystem->Exit(1);
0340   return;
0341 }