g4simulation/g4micromegas/PHG4MicromegasSteppingAction.cc

0001 /*!
0002  * \file PHG4MicromegasSteppingAction.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005
0006 #include "PHG4MicromegasSteppingAction.h"
0007
0008 #include "PHG4MicromegasDetector.h"
0009
0010 #include <phparameter/PHParameters.h>
0011
0012 #include <g4detectors/PHG4StepStatusDecode.h>
0013
0014 #include <g4main/PHG4Hit.h>
0015 #include <g4main/PHG4HitContainer.h>
0016 #include <g4main/PHG4Hitv1.h>
0017 #include <g4main/PHG4Shower.h>
0018 #include <g4main/PHG4SteppingAction.h>
0019 #include <g4main/PHG4TrackUserInfoV1.h>
0020
0021 #include <phool/getClass.h>
0022
0023 #include <TSystem.h>
0024
0025 #include <Geant4/G4ParticleDefinition.hh>
0026 #include <Geant4/G4ReferenceCountedHandle.hh>
0027 #include <Geant4/G4Step.hh>
0028 #include <Geant4/G4StepPoint.hh>
0029 #include <Geant4/G4StepStatus.hh>
0030 #include <Geant4/G4String.hh>
0031 #include <Geant4/G4SystemOfUnits.hh>
0032 #include <Geant4/G4ThreeVector.hh>
0033 #include <Geant4/G4TouchableHandle.hh>
0034 #include <Geant4/G4Track.hh>
0035 #include <Geant4/G4TrackStatus.hh>
0036 #include <Geant4/G4Types.hh>
0037 #include <Geant4/G4VPhysicalVolume.hh>
0038 #include <Geant4/G4VTouchable.hh>
0039 #include <Geant4/G4VUserTrackInformation.hh>
0040
0041 #include <cmath>
0042 #include <iostream>
0043 #include <string>
0044
0045 class PHCompositeNode;
0046
0047 //____________________________________________________________________________..
0048 PHG4MicromegasSteppingAction::PHG4MicromegasSteppingAction(PHG4MicromegasDetector *detector, const PHParameters *parameters)
0049   : PHG4SteppingAction(detector->GetName())
0050   , m_Detector(detector)
0051   , m_Params(parameters)
0052   , m_ActiveFlag(m_Params->get_int_param("active"))
0053   , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
0054 {
0055 }
0056
0057 //____________________________________________________________________________..
0058 // This is the implementation of the G4 UserSteppingAction
0059 bool PHG4MicromegasSteppingAction::UserSteppingAction(const G4Step *aStep, bool /*was_used*/)
0060 {
0061   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0062   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0063
0064   // get volume of the current step
0065   G4VPhysicalVolume *volume = touch->GetVolume();
0066   int whichactive = m_Detector->IsInDetector(volume);
0067   if (whichactive == 0)
0068   {
0069     return false;
0070   }
0071
0072   // collect energy and track length step by step
0073   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0074   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0075   const G4Track *aTrack = aStep->GetTrack();
0076
0077   // if this detector stops everything, just put all kinetic energy into edep
0078   if (m_BlackHoleFlag)
0079   {
0080     edep = aTrack->GetKineticEnergy() / GeV;
0081     auto killtrack = const_cast<G4Track *>(aTrack);
0082     killtrack->SetTrackStatus(fStopAndKill);
0083   }
0084
0085   bool geantino =
0086       aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0087       aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos;
0088
0089   G4StepPoint *prePoint = aStep->GetPreStepPoint();
0090   G4StepPoint *postPoint = aStep->GetPostStepPoint();
0091
0092   // Here we have to decide if we need to create a new hit.  Normally this should
0093   // only be neccessary if a G4 Track enters a new volume or is freshly created
0094   // For this we look at the step status of the prePoint (beginning of the G4 Step).
0095   // This should be either fGeomBoundary (G4 Track crosses into volume) or
0096   // fUndefined (G4 Track newly created)
0097   // Sadly over the years with different G4 versions we have observed cases where
0098   // G4 produces "impossible hits" which we try to catch here
0099   // These errors were always rare and it is not clear if they still exist but we
0100   // still check for them for safety. We can reproduce G4 runs identically (if given
0101   // the sequence of random number seeds you find in the log), the printouts help
0102   // us giving the G4 support information about those failures
0103   switch (prePoint->GetStepStatus())
0104   {
0105   case fPostStepDoItProc:
0106     if (m_SavePostStepStatus == fGeomBoundary)
0107     {
0108       // this is an impossible G4 Step print out diagnostic to help debug, not sure if
0109       // this is still with us
0110       std::cout << "PHG4MicromegasSteppingAction::UserSteppingAction - " << GetName() << ": New Hit for  " << std::endl;
0111       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0112                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0113                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0114                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
0115                 << std::endl;
0116
0117       std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0118       std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0119       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0120     }
0121     break;
0122
0123   // These are the normal cases
0124   case fGeomBoundary:
0125   case fUndefined:
0126   {
0127     if (!m_hit)
0128     {
0129       m_hit.reset(new PHG4Hitv1());
0130     }
0131
0132     if (whichactive > 0)
0133     {
0134       // assign layer
0135       m_hit->set_layer(m_Detector->get_layer(volume));
0136
0137       const auto tileid = m_Detector->get_tileid(volume);
0138       m_hit->set_property(PHG4Hit::prop_index_i, tileid);
0139       // momentum
0140       m_hit->set_px(0, prePoint->GetMomentum().x() / GeV);
0141       m_hit->set_py(0, prePoint->GetMomentum().y() / GeV);
0142       m_hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
0143
0144       m_hit->set_eion(0);
0145       m_SaveHitContainer = m_HitContainer;
0146     }
0147     else
0148     {
0149       m_SaveHitContainer = m_SupportHitContainer;
0150     }
0151     // here we set the entrance values in cm
0152     m_hit->set_x(0, prePoint->GetPosition().x() / cm);
0153     m_hit->set_y(0, prePoint->GetPosition().y() / cm);
0154     m_hit->set_z(0, prePoint->GetPosition().z() / cm);
0155
0156     // time in ns
0157     m_hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0158
0159     // set the track ID
0160     m_hit->set_trkid(aTrack->GetTrackID());
0161     m_SaveTrackId = aTrack->GetTrackID();
0162
0163     // reset the initial energy deposit
0164     m_EdepSum = 0;
0165     m_EionSum = 0;
0166     m_hit->set_edep(0);
0167
0168     // this is for the tracking of the truth info
0169     if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0170     {
0171       if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0172       {
0173         m_hit->set_trkid(pp->GetUserTrackId());
0174         pp->GetShower()->add_g4hit_id(m_SaveHitContainer->GetID(), m_hit->get_hit_id());
0175       }
0176     }
0177     break;
0178   }
0179
0180   default:
0181     break;
0182   }
0183
0184   if (!m_hit || !std::isfinite(m_hit->get_x(0)))
0185   {
0186     std::cout << GetName() << ": hit was not created" << std::endl;
0187     std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0188               << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0189               << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0190               << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus)
0191               << std::endl;
0192     std::cout << "last track: " << m_SaveTrackId << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0193     std::cout << "phys pre vol: " << volume->GetName() << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0194     std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName() << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0195
0196     // This is fatal - a hit from nowhere. This needs to be looked at and fixed
0197     gSystem->Exit(1);
0198   }
0199
0200   // check if track id matches the initial one when the hit was created
0201   if (aTrack->GetTrackID() != m_SaveTrackId)
0202   {
0203     std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0204     std::cout << "saved track: " << m_SaveTrackId
0205               << ", current trackid: " << aTrack->GetTrackID()
0206               << ", prestep status: " << prePoint->GetStepStatus()
0207               << ", previous post step status: " << m_SavePostStepStatus << std::endl;
0208     // This is fatal - a hit from nowhere. This needs to be looked at and fixed
0209     gSystem->Exit(1);
0210   }
0211
0212   // We need to cache a few things from one step to the next
0213   // to identify impossible hits and subsequent debugging printout
0214   m_SavePreStepStatus = prePoint->GetStepStatus();
0215   m_SavePostStepStatus = postPoint->GetStepStatus();
0216   m_SaveVolPre = volume;
0217   m_SaveVolPost = touchpost->GetVolume();
0218
0219   // here we just update the exit values, it will be overwritten
0220   // for every step until we leave the volume or the particle
0221   // ceases to exist
0222   // sum up the energy to get total deposited
0223   m_EdepSum += edep;
0224   m_EionSum += eion;
0225
0226   // if any of these conditions is true this is the last step in
0227   // this volume and we need to save the hit
0228   // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0229   // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0230   // (happens when your detector goes outside world volume)
0231   // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0232   // aTrack->GetTrackStatus() == fStopAndKill is also set)
0233   // aTrack->GetTrackStatus() == fStopAndKill: track ends
0234   if (postPoint->GetStepStatus() == fGeomBoundary ||
0235       postPoint->GetStepStatus() == fWorldBoundary ||
0236       postPoint->GetStepStatus() == fAtRestDoItProc ||
0237       aTrack->GetTrackStatus() == fStopAndKill)
0238   {
0239     // save only hits with energy deposit (or geantino)
0240     if (m_EdepSum > 0 || geantino)
0241     {
0242       // update values at exit coordinates and set keep flag
0243       // of track to keep
0244       m_hit->set_x(1, postPoint->GetPosition().x() / cm);
0245       m_hit->set_y(1, postPoint->GetPosition().y() / cm);
0246       m_hit->set_z(1, postPoint->GetPosition().z() / cm);
0247
0248       m_hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0249       if (whichactive > 0)
0250       {
0251         m_hit->set_px(1, postPoint->GetMomentum().x() / GeV);
0252         m_hit->set_py(1, postPoint->GetMomentum().y() / GeV);
0253         m_hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
0254       }
0255
0256       if (G4VUserTrackInformation *p = aTrack->GetUserInformation())
0257       {
0258         if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
0259         {
0260           pp->SetKeep(1);  // we want to keep the track
0261         }
0262       }
0263
0264       if (geantino)
0265       {
0266         m_hit->set_edep(-1);
0267       }
0268       else
0269       {
0270         m_hit->set_edep(m_EdepSum);
0271       }
0272       if (whichactive > 0)
0273       {
0274         if (geantino)
0275         {
0276           m_hit->set_eion(-1);
0277         }
0278         else
0279         {
0280           m_hit->set_eion(m_EionSum);
0281         }
0282       }
0283
0284       // add in container
0285       m_SaveHitContainer->AddHit(m_hit->get_layer(), m_hit.get());
0286
0287       // ownership is transferred to container
0288       // so we release the hit
0289       // cppcheck says to check the return code, likely it should be non-zero
0290       // cppcheck-suppress *
0291       m_hit.release();
0292     }
0293     else
0294     {
0295       m_hit->Reset();
0296     }
0297   }
0298   // return true to indicate the hit was used
0299   return true;
0300 }
0301
0302 //____________________________________________________________________________..
0303 void PHG4MicromegasSteppingAction::SetInterfacePointers(PHCompositeNode *topNode)
0304 {
0305   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0306   m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
0307
0308   if (!m_HitContainer)
0309   {
0310     if (!m_Params->get_int_param("blackhole"))  // not messed up if we have a black hole
0311     {
0312       std::cout << "PHG4MicromegasSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0313       gSystem->Exit(1);
0314     }
0315   }
0316   // this is perfectly fine if support hits are disabled
0317   if (!m_SupportHitContainer)
0318   {
0319     if (Verbosity() > 0)
0320     {
0321       std::cout << "PHG4MicromegasSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
0322     }
0323   }
0324 }
0325
0326 void PHG4MicromegasSteppingAction::SetHitNodeName(const std::string &type, const std::string &name)
0327 {
0328   if (type == "G4HIT")
0329   {
0330     m_HitNodeName = name;
0331     return;
0332   }
0333   else if (type == "G4HIT_SUPPORT")
0334   {
0335     m_SupportNodeName = name;
0336     return;
0337   }
0338   std::cout << "Invalid output hit node type " << type << std::endl;
0339   gSystem->Exit(1);
0340   return;
0341 }