g4simulation/g4tpc/PHG4TpcSteppingAction.cc

0001 #include "PHG4TpcSteppingAction.h"
0002 #include "PHG4TpcDetector.h"
0003
0004 #include <g4detectors/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
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>  // for G4ReferenceCountedHandle
0022 #include <Geant4/G4Step.hh>
0023 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0024 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fPostSt...
0025 #include <Geant4/G4String.hh>      // for G4String
0026 #include <Geant4/G4SystemOfUnits.hh>
0027 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0028 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0029 #include <Geant4/G4Track.hh>                  // for G4Track
0030 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0031 #include <Geant4/G4Types.hh>                  // for G4double
0032 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0033 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0034 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0035
0036 #include <cmath>    // for isfinite
0037 #include <cstdlib>  // for exit
0038 #include <iostream>
0039 #include <string>  // for operator<<, operator+
0040
0041 class PHCompositeNode;
0042
0043 //____________________________________________________________________________..
0044 PHG4TpcSteppingAction::PHG4TpcSteppingAction(PHG4TpcDetector* detector, const PHParameters* parameters)
0045   : PHG4SteppingAction(detector->GetName())
0046   , m_Detector(detector)
0047   , m_Params(parameters)
0048   , m_IsBlackHoleFlag(m_Params->get_int_param("blackhole"))
0049 {
0050   if (std::isfinite(m_Params->get_double_param("steplimits")))
0051   {
0052     m_UseG4StepsFlag = 1;
0053   }
0054   SetName(m_Detector->GetName());
0055 }
0056
0057 PHG4TpcSteppingAction::~PHG4TpcSteppingAction()
0058 {
0059   // if the last hit was a zero energie deposit hit, it is just reset
0060   // and the memory is still allocated, so we need to delete it here
0061   // if the last hit was saved, hit is a nullptr pointer which are
0062   // legal to delete (it results in a no operation)
0063   delete m_Hit;
0064 }
0065 //____________________________________________________________________________..
0066 bool PHG4TpcSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0067 {
0068   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0069   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0070   // get volume of the current step
0071   G4VPhysicalVolume* volume = touch->GetVolume();
0072
0073   // m_Detector->IsInTpc(volume)
0074   // returns
0075   //  0 is outside of Tpc
0076   //  1 is inside tpc gas
0077   // <0 is in tpc support structure (cage, endcaps,...)
0078
0079   int whichactive = m_Detector->IsInTpc(volume);
0080   if (!whichactive)
0081   {
0082     return false;
0083   }
0084   // collect energy and track length step by step
0085   G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0086   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0087   const G4Track* aTrack = aStep->GetTrack();
0088
0089   // if this block stops everything, just put all kinetic energy into edep
0090   if (m_IsBlackHoleFlag)
0091   {
0092     edep = aTrack->GetKineticEnergy() / GeV;
0093     G4Track* killtrack = const_cast<G4Track*>(aTrack);
0094     killtrack->SetTrackStatus(fStopAndKill);
0095   }
0096
0097   bool geantino = false;
0098
0099   // the check for the pdg code speeds things up, I do not want to make
0100   // an expensive string compare for every track when we know
0101   // geantino or chargedgeantino has pid=0
0102   if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0103       aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0104   {
0105     geantino = true;
0106   }
0107   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0108   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0109   int prepointstatus = prePoint->GetStepStatus();
0110
0111   //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0112   //       std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0113   //       std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0114
0115   if ((m_UseG4StepsFlag > 0 && whichactive > 0) ||
0116       prepointstatus == fGeomBoundary ||
0117       prepointstatus == fUndefined ||
0118       (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary))
0119   {
0120     unsigned int layer_id = 99;  // no layer number for the hit, use a non-existent one for now, replace it later
0121     // this is for debugging weird occurances we have occasionally
0122     if (prepointstatus == fPostStepDoItProc && m_SavePostStepStatus == fGeomBoundary)
0123     {
0124       std::cout << GetName() << ": New Hit for  " << std::endl;
0125       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0126                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0127                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0128                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0129       std::cout << "last track: " << m_SaveTrackId
0130                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0131       std::cout << "phys pre vol: " << volume->GetName()
0132                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0133       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0134                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0135     }
0136     // if previous hit was saved, hit pointer was set to nullptr
0137     // and we have to make a new one
0138     if (!m_Hit)
0139     {
0140       m_Hit = new PHG4Hitv1();
0141     }
0142     m_Hit->set_layer(layer_id);
0143     // here we set the entrance values in cm
0144     m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
0145     m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
0146     m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
0147
0148     // momentum
0149     m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
0150     m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
0151     m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
0152
0153     // time in ns
0154     m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0155     // set and save the track ID
0156     m_Hit->set_trkid(aTrack->GetTrackID());
0157     m_SaveTrackId = aTrack->GetTrackID();
0158     // set the initial energy deposit
0159     m_Hit->set_edep(0);
0160     if (whichactive > 0)  // return of IsInTpcDetector, > 0 hit in tpc gas volume, < 0 hit in support structures
0161     {
0162       m_Hit->set_eion(0);
0163       // Now save the container we want to add this hit to
0164       m_CurrentHitContainer = m_HitContainer;
0165     }
0166     else
0167     {
0168       m_CurrentHitContainer = m_AbsorberHitContainer;
0169     }
0170     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0171     {
0172       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0173       {
0174         m_Hit->set_trkid(pp->GetUserTrackId());
0175         m_Hit->set_shower_id(pp->GetShower()->get_id());
0176         m_Shower = pp->GetShower();
0177       }
0178     }
0179     // some sanity checks for inconsistencies
0180     // check if this hit was created, if not print out last post step status
0181     if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
0182     {
0183       std::cout << GetName() << ": hit was not created" << std::endl;
0184       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0185                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0186                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePreStepStatus)
0187                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(m_SavePostStepStatus) << std::endl;
0188       std::cout << "last track: " << m_SaveTrackId
0189                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0190       std::cout << "phys pre vol: " << volume->GetName()
0191                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0192       std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
0193                 << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
0194       exit(1);
0195     }
0196     // check if track id matches the initial one when the hit was created
0197     if (aTrack->GetTrackID() != m_SaveTrackId)
0198     {
0199       std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0200       std::cout << "saved track: " << m_SaveTrackId
0201                 << ", current trackid: " << aTrack->GetTrackID()
0202                 << ", prestep status: " << prePoint->GetStepStatus()
0203                 << ", previous post step status: " << m_SavePostStepStatus
0204                 << std::endl;
0205
0206       exit(1);
0207     }
0208     m_SavePreStepStatus = prePoint->GetStepStatus();
0209     m_SavePostStepStatus = postPoint->GetStepStatus();
0210     m_SaveVolPre = volume;
0211     m_SaveVolPost = touchpost->GetVolume();
0212     // here we just update the exit values, it will be overwritten
0213     // for every step until we leave the volume or the particle
0214     // ceases to exist
0215     m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
0216     m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
0217     m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
0218
0219     m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
0220     m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
0221     m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
0222
0223     m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0224
0225     // sum up the energy to get total deposited
0226     m_Hit->set_edep(m_Hit->get_edep() + edep);
0227     if (whichactive > 0)
0228     {
0229       m_Hit->set_eion(m_Hit->get_eion() + eion);
0230     }
0231     if (geantino)
0232     {
0233       m_Hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0234       if (whichactive > 0)
0235       {
0236         m_Hit->set_eion(-1);
0237       }
0238     }
0239     if (edep > 0)
0240     {
0241       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0242       {
0243         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0244         {
0245           pp->SetKeep(1);  // we want to keep the track
0246         }
0247       }
0248     }
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 ((m_UseG4StepsFlag > 0 && whichactive > 0) ||
0259         postPoint->GetStepStatus() == fGeomBoundary ||
0260         postPoint->GetStepStatus() == fWorldBoundary ||
0261         postPoint->GetStepStatus() == fAtRestDoItProc ||
0262         aTrack->GetTrackStatus() == fStopAndKill)
0263     {
0264       // save only hits with energy deposit (or -1 for geantino)
0265       if (m_Hit->get_edep())
0266       {
0267         m_CurrentHitContainer->AddHit(layer_id, m_Hit);
0268         if (m_Shower)
0269         {
0270           m_Shower->add_g4hit_id(m_CurrentHitContainer->GetID(), m_Hit->get_hit_id());
0271         }
0272         // promote to double to force double sqrt
0273         double rin = sqrt((double) (m_Hit->get_x(0) * m_Hit->get_x(0) + m_Hit->get_y(0) * m_Hit->get_y(0)));
0274         double rout = sqrt((double) (m_Hit->get_x(1) * m_Hit->get_x(1) + m_Hit->get_y(1) * m_Hit->get_y(1)));
0275         if (Verbosity() > 10)
0276         {
0277           if ((rin > 69.0 && rin < 70.125) || (rout > 69.0 && rout < 70.125))
0278           {
0279             std::cout << "Added Tpc g4hit with rin, rout = " << rin << "  " << rout
0280                       << " g4hitid " << m_Hit->get_hit_id() << std::endl;
0281             std::cout << " xin " << m_Hit->get_x(0)
0282                       << " yin " << m_Hit->get_y(0)
0283                       << " zin " << m_Hit->get_z(0)
0284                       << " rin " << rin
0285                       << std::endl;
0286             std::cout << " xout " << m_Hit->get_x(1)
0287                       << " yout " << m_Hit->get_y(1)
0288                       << " zout " << m_Hit->get_z(1)
0289                       << " rout " << rout
0290                       << std::endl;
0291             std::cout << " xav " << (m_Hit->get_x(1) + m_Hit->get_x(0)) / 2.0
0292                       << " yav " << (m_Hit->get_y(1) + m_Hit->get_y(0)) / 2.0
0293                       << " zav " << (m_Hit->get_z(1) + m_Hit->get_z(0)) / 2.0
0294                       << " rav " << (rout + rin) / 2.0
0295                       << std::endl;
0296           }
0297         }
0298
0299         // ownership has been transferred to container, set to null
0300         // so we will create a new hit for the next track
0301         m_Hit = nullptr;
0302       }
0303       else
0304       {
0305         // if this hit has no energy deposit, just reset it for reuse
0306         // this means we have to delete it in the dtor. If this was
0307         // the last hit we processed the memory is still allocated
0308         m_Hit->Reset();
0309       }
0310     }
0311     // return true to indicate the hit was used
0312     return true;
0313   }
0314   else
0315   {
0316     return false;
0317   }
0318 }
0319
0320 //____________________________________________________________________________..
0321 void PHG4TpcSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0322 {
0323   m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
0324   m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
0325
0326   // if we do not find the node it's messed up.
0327   if (!m_HitContainer)
0328   {
0329     std::cout << "PHG4TpcSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
0330   }
0331   if (!m_AbsorberHitContainer)
0332   {
0333     if (Verbosity() > 1)
0334     {
0335       std::cout << "PHG4HcalSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
0336     }
0337   }
0338 }
0339
0340 void PHG4TpcSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
0341 {
0342   if (type == "G4HIT")
0343   {
0344     m_HitNodeName = name;
0345     return;
0346   }
0347   else if (type == "G4HIT_ABSORBER")
0348   {
0349     m_AbsorberNodeName = name;
0350     return;
0351   }
0352   std::cout << "Invalid output hit node type " << type << std::endl;
0353   gSystem->Exit(1);
0354   return;
0355 }