g4simulation/g4detectors/BeamLineMagnetSteppingAction.cc

0001 #include "BeamLineMagnetSteppingAction.h"
0002
0003 #include "BeamLineMagnetDetector.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 <TSystem.h>
0018
0019 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0020 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0021 #include <Geant4/G4Step.hh>
0022 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0023 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRest...
0024 #include <Geant4/G4String.hh>      // for G4String
0025 #include <Geant4/G4SystemOfUnits.hh>
0026 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0027 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0028 #include <Geant4/G4Track.hh>                  // for G4Track
0029 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0030 #include <Geant4/G4Types.hh>                  // for G4double
0031 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0032 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0033 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0034
0035 #include <cmath>  // for isfinite
0036 #include <iostream>
0037 #include <string>  // for operator<<, string
0038
0039 class PHCompositeNode;
0040
0041 //____________________________________________________________________________..
0042 BeamLineMagnetSteppingAction::BeamLineMagnetSteppingAction(BeamLineMagnetDetector* detector, const PHParameters* parameters)
0043   : PHG4SteppingAction(detector->GetName())
0044   , m_Detector(detector)
0045   , m_Params(parameters)
0046   , m_ActiveFlag(m_Params->get_int_param("active"))
0047   , m_AbsorberActiveFlag(m_Params->get_int_param("absorberactive"))
0048   , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
0049 {
0050 }
0051
0052 BeamLineMagnetSteppingAction::~BeamLineMagnetSteppingAction()
0053 {
0054   // if the last hit was a zero energie deposit hit, it is just reset
0055   // and the memory is still allocated, so we need to delete it here
0056   // if the last hit was saved, hit is a nullptr pointer which are
0057   // legal to delete (it results in a no operation)
0058   delete m_Hit;
0059 }
0060
0061 //____________________________________________________________________________..
0062 bool BeamLineMagnetSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0063 {
0064   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0065   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0066   // get volume of the current step
0067   G4VPhysicalVolume* volume = touch->GetVolume();
0068   // IsInVTX(volume) returns
0069   //  == 0 outside of pstof
0070   //   > 0 for hits in active volume
0071   //  < 0 for hits in passive material
0072   int whichactive = m_Detector->IsInBeamLineMagnet(volume);
0073   if (!whichactive)
0074   {
0075     return false;
0076   }
0077
0078   // collect energy and track length step by step
0079   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0080   const G4Track* aTrack = aStep->GetTrack();
0081
0082   // if this block stops everything ONLY if in the magnet iron, just put all kinetic energy into edep
0083   if (whichactive == -1)
0084   {
0085     if (m_BlackHoleFlag)
0086     {
0087       edep = aTrack->GetKineticEnergy() / GeV;
0088       G4Track* killtrack = const_cast<G4Track*>(aTrack);
0089       killtrack->SetTrackStatus(fStopAndKill);
0090     }
0091     if (!m_AbsorberActiveFlag)
0092     {
0093       return false;
0094     }
0095   }
0096   if (whichactive == -2)
0097   {
0098     edep = aTrack->GetKineticEnergy() / GeV;
0099     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0100     killtrack->SetTrackStatus(fStopAndKill);
0101   }
0102   if (!m_ActiveFlag)
0103   {
0104     return false;
0105   }
0106
0107   int magnet_id = volume->GetCopyNo();  // magnet id is stored in copy number
0108   bool geantino = false;
0109   // the check for the pdg code speeds things up, I do not want to make
0110   // an expensive string compare for every track when we know
0111   // geantino or chargedgeantino has pid=0
0112   if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0113       aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0114   {
0115     geantino = true;
0116   }
0117   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0118   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0119   //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0120   //       std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0121   //       std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0122
0123   switch (prePoint->GetStepStatus())
0124   {
0125   case fPostStepDoItProc:
0126     if (m_SavePostStepStatus != fGeomBoundary)
0127     {
0128       break;
0129     }
0130     else
0131     {
0132       std::cout << GetName() << ": New Hit for  " << std::endl;
0133       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0134                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0135                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0136                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0137       std::cout << "last track: " << m_SaveTrackId
0138                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0139       std::cout << "phys pre vol: " << volume->GetName()
0140                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0141       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0142                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0143     }
0144     [[fallthrough]];
0145   case fGeomBoundary:
0146   case fUndefined:
0147     if (!m_Hit)
0148     {
0149       m_Hit = new PHG4Hitv1();
0150     }
0151     m_Hit->set_layer(magnet_id);
0152     // here we set the entrance values in cm
0153     m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0154     m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0155     m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0156     // time in ns
0157     m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0158     // set the track ID
0159     m_Hit->set_trkid(aTrack->GetTrackID());
0160     m_SaveTrackId = aTrack->GetTrackID();
0161     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0162     {
0163       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0164       {
0165         m_Hit->set_trkid(pp->GetUserTrackId());
0166       }
0167     }
0168     // set the initial energy deposit
0169     m_EdepSum = 0;
0170     if (whichactive > 0)
0171     {
0172       m_SaveHitContainer = m_HitContainer;
0173     }
0174     else
0175     {
0176       m_SaveHitContainer = m_AbsorberHitContainer;
0177     }
0178     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0179     {
0180       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0181       {
0182         m_Hit->set_trkid(pp->GetUserTrackId());
0183         pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_Hit->get_hit_id());
0184       }
0185     }
0186
0187     break;
0188   default:
0189     break;
0190   }
0191
0192   // some sanity checks for inconsistencies
0193   // check if this hit was created, if not print out last post step status
0194   if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0195   {
0196     std::cout << GetName() << ": hit was not created" << std::endl;
0197     std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0198               << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0199               << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0200               << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0201     std::cout << "last track: " << m_SaveTrackId
0202               << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0203     std::cout << "phys pre vol: " << volume->GetName()
0204               << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0205     std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0206               << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0207     gSystem->Exit(1);
0208   }
0209   // check if track id matches the initial one when the hit was created
0210   if (aTrack->GetTrackID() != m_SaveTrackId)
0211   {
0212     std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0213     std::cout << "saved track: " << m_SaveTrackId
0214               << ", current trackid: " << aTrack->GetTrackID()
0215               << ", prestep status: " << prePoint->GetStepStatus()
0216               << ", previous post step status: " << m_SavePostStepStatus
0217               << std::endl;
0218
0219     gSystem->Exit(1);
0220   }
0221   m_SavePreStepStatus = prePoint->GetStepStatus();
0222   m_SavePostStepStatus = postPoint->GetStepStatus();
0223   m_SaveVolPre = volume;
0224   m_SaveVolPost = touchpost->GetVolume();
0225
0226   // here we just update the exit values, it will be overwritten
0227   // for every step until we leave the volume or the particle
0228   // ceases to exist
0229   // sum up the energy to get total deposited
0230   m_EdepSum += edep;
0231   // if any of these conditions is true this is the last step in
0232   // this volume and we need to save the hit
0233   // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0234   // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0235   // (happens when your detector goes outside world volume)
0236   // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0237   // aTrack->GetTrackStatus() == fStopAndKill is also set)
0238   // aTrack->GetTrackStatus() == fStopAndKill: track ends
0239   if (postPoint->GetStepStatus() == fGeomBoundary ||
0240       postPoint->GetStepStatus() == fWorldBoundary ||
0241       postPoint->GetStepStatus() == fAtRestDoItProc ||
0242       aTrack->GetTrackStatus() == fStopAndKill)
0243   {
0244     // save only hits with energy deposit (or geantino)
0245     if (m_EdepSum > 0 || geantino)
0246     {
0247       // update values at exit coordinates and set keep flag
0248       // of track to keep
0249       m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0250       m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0251       m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0252
0253       m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0254       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0255       {
0256         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0257         {
0258           pp->SetKeep(1);  // we want to keep the track
0259         }
0260       }
0261       if (geantino)
0262       {
0263         m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0264       }
0265       else
0266       {
0267         m_Hit->set_edep(m_EdepSum);
0268       }
0269       m_SaveHitContainer->AddHit(magnet_id, m_Hit);
0270       // ownership has been transferred to container, set to null
0271       // so we will create a new hit for the next track
0272       m_Hit = nullptr;
0273     }
0274     else
0275     {
0276       // if this hit has no energy deposit, just reset it for reuse
0277       // this means we have to delete it in the dtor. If this was
0278       // the last hit we processed the memory is still allocated
0279       m_Hit->Reset();
0280     }
0281   }
0282   // return true to indicate the hit was used
0283   return true;
0284 }
0285
0286 //____________________________________________________________________________..
0287 void BeamLineMagnetSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0288 {
0289   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0290   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
0291
0292   // If this is a black hole but otherwise not an active detector there is
0293   // no hitnode, otherwise this will crash in the user stepping action
0294   if (!m_HitContainer)
0295   {
0296     if (!m_BlackHoleFlag)
0297     {
0298       std::cout << "BeamLineMagnetSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0299     }
0300   }
0301   if (!m_AbsorberHitContainer)
0302   {
0303     if (Verbosity() > 1)
0304     {
0305       std::cout << "BeamLineMagnetSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
0306     }
0307   }
0308 }
0309
0310 void BeamLineMagnetSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
0311 {
0312   if (type == "G4HIT")
0313   {
0314     m_HitNodeName = name;
0315     return;
0316   }
0317   else if (type == "G4HIT_ABSORBER")
0318   {
0319     m_AbsorberNodeName = name;
0320     return;
0321   }
0322   std::cout << "Invalid output hit node type " << type << std::endl;
0323   gSystem->Exit(1);
0324   return;
0325 }