g4simulation/g4detectors/PHG4CylinderSteppingAction.cc

0001 #include "PHG4CylinderSteppingAction.h"
0002 #include "PHG4CylinderDetector.h"
0003 #include "PHG4CylinderSubsystem.h"
0004
0005 #include "PHG4StepStatusDecode.h"
0006
0007 #include <phparameter/PHParameters.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>  // for PHG4SteppingAction
0014
0015 #include <g4main/PHG4TrackUserInfoV1.h>
0016
0017 #include <phool/getClass.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, fPostSt...
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 #pragma GCC diagnostic push
0036 #pragma GCC diagnostic ignored "-Wshadow"
0037 #include <boost/io/ios_state.hpp>
0038 #pragma GCC diagnostic pop
0039
0040 #include <cmath>    // for isfinite, copysign
0041 #include <cstdlib>  // for exit
0042 #include <iomanip>
0043 #include <iostream>
0044 #include <string>  // for operator<<, char_traits
0045
0046 class PHCompositeNode;
0047
0048 //____________________________________________________________________________..
0049 PHG4CylinderSteppingAction::PHG4CylinderSteppingAction(PHG4CylinderSubsystem* subsys, PHG4CylinderDetector* detector, const PHParameters* parameters)
0050   : PHG4SteppingAction(detector->GetName())
0051   , m_Subsystem(subsys)
0052   , m_Detector(detector)
0053   , m_Params(parameters)
0054   , m_HitContainer(nullptr)
0055   , m_Hit(nullptr)
0056   , m_SaveShower(nullptr)
0057   , m_SaveVolPre(nullptr)
0058   , m_SaveVolPost(nullptr)
0059   , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
0060   , m_SaveTrackId(-1)
0061   , m_SavePreStepStatus(-1)
0062   , m_SavePostStepStatus(-1)
0063   , m_ActiveFlag(m_Params->get_int_param("active"))
0064   , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
0065   , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
0066   , m_Zmin(m_Params->get_double_param("place_z") * cm - m_Params->get_double_param("length") * cm / 2.)
0067   , m_Zmax(m_Params->get_double_param("place_z") * cm + m_Params->get_double_param("length") * cm / 2.)
0068   , m_Tmin(m_Params->get_double_param("tmin") * ns)
0069   , m_Tmax(m_Params->get_double_param("tmax") * ns)
0070 {
0071   // G4 seems to have issues in the um range
0072   m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
0073   m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
0074 }
0075
0076 PHG4CylinderSteppingAction::~PHG4CylinderSteppingAction()
0077 {
0078   // if the last hit was a zero energie deposit hit, it is just reset
0079   // and the memory is still allocated, so we need to delete it here
0080   // if the last hit was saved, hit is a nullptr pointer which are
0081   // legal to delete (it results in a no operation)
0082   delete m_Hit;
0083 }
0084
0085 //____________________________________________________________________________..
0086 bool PHG4CylinderSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0087 {
0088   // get volume of the current step
0089   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0090   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0091
0092   G4VPhysicalVolume* volume = touch->GetVolume();
0093   // G4 just calls  UserSteppingAction for every step (and we loop then over all our
0094   // steppingactions. First we have to check if we are actually in our volume
0095   if (!m_Detector->IsInCylinder(volume))
0096   {
0097     return false;
0098   }
0099
0100   // collect energy and track length step by step
0101   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0102
0103   const G4Track* aTrack = aStep->GetTrack();
0104
0105   // if this cylinder stops everything, just put all kinetic energy into edep
0106   if (m_BlackHoleFlag)
0107   {
0108     if ((!std::isfinite(m_Tmin) && !std::isfinite(m_Tmax)) ||
0109         aTrack->GetGlobalTime() < m_Tmin ||
0110         aTrack->GetGlobalTime() > m_Tmax)
0111     {
0112       edep = aTrack->GetKineticEnergy() / GeV;
0113       G4Track* killtrack = const_cast<G4Track*>(aTrack);
0114       killtrack->SetTrackStatus(fStopAndKill);
0115     }
0116   }
0117
0118   int layer_id = m_Detector->get_Layer();
0119   // test if we are active
0120   if (m_ActiveFlag)
0121   {
0122     bool geantino = false;
0123     // the check for the pdg code speeds things up, I do not want to make
0124     // an expensive string compare for every track when we know
0125     // geantino or chargedgeantino has pid=0
0126     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0127         aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0128     {
0129       geantino = true;
0130     }
0131     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0132     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0133     //        std::cout << "time prepoint: " << prePoint->GetGlobalTime()/ns << std::endl;
0134     //        std::cout << "time postpoint: " << postPoint->GetGlobalTime()/ns << std::endl;
0135     //        std::cout << "kinetic energy: " <<  aTrack->GetKineticEnergy()/GeV << std::endl;
0136     //       G4ParticleDefinition* def = aTrack->GetDefinition();
0137     //       std::cout << "Particle: " << def->GetParticleName() << std::endl;
0138     int prepointstatus = prePoint->GetStepStatus();
0139     if (prepointstatus == fGeomBoundary ||
0140         prepointstatus == fUndefined ||
0141         (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary) ||
0142         m_UseG4StepsFlag > 0)
0143     {
0144       if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
0145       {
0146         std::cout << GetName() << ": New Hit for  " << std::endl;
0147         std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0148                   << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0149                   << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0150                   << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0151         std::cout << "last track: " << m_SaveTrackId
0152                   << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0153         std::cout << "phys pre vol: " << volume->GetName()
0154                   << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0155         std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0156                   << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0157       }
0158
0159       if (!m_Hit)
0160       {
0161         m_Hit = new PHG4Hitv1();
0162       }
0163
0164       m_Hit->set_layer((unsigned int) layer_id);
0165
0166       // here we set the entrance values in cm
0167       m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0168       m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0169       m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0170
0171       m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
0172       m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
0173       m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
0174
0175       // time in ns
0176       m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0177       // set and save the track ID
0178       m_Hit->set_trkid(aTrack->GetTrackID());
0179       m_SaveTrackId = aTrack->GetTrackID();
0180       // set the initial energy deposit
0181       m_Hit->set_edep(0);
0182       if (!geantino && !m_BlackHoleFlag)
0183       {
0184         m_Hit->set_eion(0);
0185       }
0186       if (m_SaveLightYieldFlag)
0187       {
0188         m_Hit->set_light_yield(0);
0189       }
0190       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0191       {
0192         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0193         {
0194           m_Hit->set_trkid(pp->GetUserTrackId());
0195           m_Hit->set_shower_id(pp->GetShower()->get_id());
0196           m_SaveShower = pp->GetShower();
0197         }
0198       }
0199
0200       if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
0201       {
0202         boost::io::ios_precision_saver ips(std::cout);
0203         std::cout << m_Detector->SuperDetector() << std::setprecision(9)
0204                   << " PHG4CylinderSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
0205                   << " outside acceptance,  zmin " << m_Zmin
0206                   << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
0207       }
0208     }
0209     // here we just update the exit values, it will be overwritten
0210     // for every step until we leave the volume or the particle
0211     // ceases to exist
0212     // some sanity checks for inconsistencies
0213     // check if this hit was created, if not print out last post step status
0214     if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0215     {
0216       std::cout << GetName() << ": hit was not created" << std::endl;
0217       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0218                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0219                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0220                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0221       std::cout << "last track: " << m_SaveTrackId
0222                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0223       std::cout << "phys pre vol: " << volume->GetName()
0224                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0225       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0226                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0227       exit(1);
0228     }
0229     m_SavePostStepStatus = postPoint->GetStepStatus();
0230     // check if track id matches the initial one when the hit was created
0231     if (aTrack->GetTrackID() != m_SaveTrackId)
0232     {
0233       std::cout << "hits do not belong to the same track" << std::endl;
0234       std::cout << "saved track: " << m_SaveTrackId
0235                 << ", current trackid: " << aTrack->GetTrackID()
0236                 << std::endl;
0237       exit(1);
0238     }
0239     m_SavePreStepStatus = prePoint->GetStepStatus();
0240     m_SavePostStepStatus = postPoint->GetStepStatus();
0241     m_SaveVolPre = volume;
0242     m_SaveVolPost = touchpost->GetVolume();
0243
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_px(1, postPoint->GetMomentum().x() / GeV);
0249     m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
0250     m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
0251
0252     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0253     // sum up the energy to get total deposited
0254     m_Hit->set_edep(m_Hit->get_edep() + edep);
0255     if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
0256     {
0257       std::cout << m_Detector->SuperDetector() << std::setprecision(9)
0258                 << " PHG4CylinderSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
0259                 << " outside acceptance zmin " << m_Zmin
0260                 << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
0261     }
0262     if (geantino)
0263     {
0264       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0265     }
0266     else
0267     {
0268       if (!m_BlackHoleFlag)
0269       {
0270         double eion = edep - aStep->GetNonIonizingEnergyDeposit() / GeV;
0271         m_Hit->set_eion(m_Hit->get_eion() + eion);
0272       }
0273     }
0274     if (m_SaveLightYieldFlag)
0275     {
0276       double light_yield = GetVisibleEnergyDeposition(aStep);
0277       m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
0278     }
0279     if (edep > 0 || m_SaveAllHitsFlag)
0280     {
0281       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0282       {
0283         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0284         {
0285           pp->SetKeep(1);  // we want to keep the track
0286         }
0287       }
0288     }
0289     // if any of these conditions is true this is the last step in
0290     // this volume and we need to save the hit
0291     // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0292     // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0293     // (happens when your detector goes outside world volume)
0294     // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0295     // aTrack->GetTrackStatus() == fStopAndKill is also set)
0296     // aTrack->GetTrackStatus() == fStopAndKill: track ends
0297     if (postPoint->GetStepStatus() == fGeomBoundary ||
0298         postPoint->GetStepStatus() == fWorldBoundary ||
0299         postPoint->GetStepStatus() == fAtRestDoItProc ||
0300         aTrack->GetTrackStatus() == fStopAndKill ||
0301         m_UseG4StepsFlag > 0)
0302     {
0303       // save only hits with energy deposit (or -1 for geantino) or if save all hits flag is set
0304       if (m_Hit->get_edep() || m_SaveAllHitsFlag)
0305       {
0306         m_HitContainer->AddHit(layer_id, m_Hit);
0307         if (m_SaveShower)
0308         {
0309           m_SaveShower->add_g4hit_id(m_HitContainer->GetID(), m_Hit->get_hit_id());
0310         }
0311         // ownership has been transferred to container, set to null
0312         // so we will create a new hit for the next track
0313         m_Hit = nullptr;
0314       }
0315       else
0316       {
0317         // if this hit has no energy deposit, just reset it for reuse
0318         // this means we have to delete it in the dtor. If this was
0319         // the last hit we processed the memory is still allocated
0320         m_Hit->Reset();
0321       }
0322     }
0323     // return true to indicate the hit was used
0324     return true;
0325   }
0326   else
0327   {
0328     return false;
0329   }
0330 }
0331
0332 //____________________________________________________________________________..
0333 void PHG4CylinderSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0334 {
0335   // Node Name is passed down from PHG4CylinderSubsystem
0336   // now look for the map and grab a pointer to it.
0337   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0338
0339   // if we do not find the node we need to scream.
0340   if (!m_HitContainer && !m_BlackHoleFlag)
0341   {
0342     std::cout << "PHG4CylinderSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0343   }
0344 }
0345
0346 bool PHG4CylinderSteppingAction::hasMotherSubsystem() const
0347 {
0348   if (m_Subsystem->GetMotherSubsystem())
0349   {
0350     return true;
0351   }
0352   return false;
0353 }