g4simulation/g4detectors/PHG4PSTOFSteppingAction.cc

0001 #include "PHG4PSTOFSteppingAction.h"
0002 #include "PHG4PSTOFDetector.h"
0003 #include "PHG4StepStatusDecode.h"
0004
0005 #include <g4main/PHG4Hit.h>
0006 #include <g4main/PHG4HitContainer.h>
0007 #include <g4main/PHG4Hitv1.h>
0008 #include <g4main/PHG4Shower.h>
0009 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0010 #include <g4main/PHG4TrackUserInfoV1.h>
0011
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>  // for PHWHERE
0014
0015 #include <TSystem.h>
0016
0017 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0018 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0019 #include <Geant4/G4Step.hh>
0020 #include <Geant4/G4StepPoint.hh>   // for G4StepPoint
0021 #include <Geant4/G4StepStatus.hh>  // for fGeomBoundary, fAtRest...
0022 #include <Geant4/G4String.hh>      // for G4String
0023 #include <Geant4/G4SystemOfUnits.hh>
0024 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0025 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0026 #include <Geant4/G4Track.hh>                  // for G4Track
0027 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0028 #include <Geant4/G4Types.hh>                  // for G4double
0029 #include <Geant4/G4VPhysicalVolume.hh>        // for G4VPhysicalVolume
0030 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0031 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0032
0033 #include <cmath>  // for isfinite
0034 #include <iostream>
0035 #include <string>  // for operator<<, string
0036
0037 class PHCompositeNode;
0038
0039 //____________________________________________________________________________..
0040 PHG4PSTOFSteppingAction::PHG4PSTOFSteppingAction(PHG4PSTOFDetector* detector, const PHParametersContainer* /*parameters*/)
0041   : PHG4SteppingAction(detector->GetName())
0042   , detector_(detector)
0043 {
0044 }
0045
0046 PHG4PSTOFSteppingAction::~PHG4PSTOFSteppingAction()
0047 {
0048   // if the last hit was a zero energie deposit hit, it is just reset
0049   // and the memory is still allocated, so we need to delete it here
0050   // if the last hit was saved, hit is a nullptr pointer which are
0051   // legal to delete (it results in a no operation)
0052   delete hit;
0053 }
0054
0055 //____________________________________________________________________________..
0056 bool PHG4PSTOFSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0057 {
0058   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0059   G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
0060   // get volume of the current step
0061   G4VPhysicalVolume* volume = touch->GetVolume();
0062   // IsInPSTOF(volume) returns
0063   //  == 0 outside of pstof
0064   //   > 0 for hits in active volume
0065   //  < 0 for hits in passive material
0066   int whichactive = detector_->IsInPSTOF(volume);
0067   if (!whichactive)
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   int layer_id = 0;  // what the heck is this?
0078   bool geantino = false;
0079   // the check for the pdg code speeds things up, I do not want to make
0080   // an expensive string compare for every track when we know
0081   // geantino or chargedgeantino has pid=0
0082   if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
0083       aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0084   {
0085     geantino = true;
0086   }
0087   G4StepPoint* prePoint = aStep->GetPreStepPoint();
0088   G4StepPoint* postPoint = aStep->GetPostStepPoint();
0089   //       std::cout << "track id " << aTrack->GetTrackID() << std::endl;
0090   //       std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
0091   //       std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
0092
0093   layer_id = touch->GetCopyNumber();
0094   if (layer_id != whichactive)
0095   {
0096     std::cout << PHWHERE << " inconsistency between G4 copy number: "
0097               << layer_id << " and module id from detector: "
0098               << whichactive << std::endl;
0099     gSystem->Exit(1);
0100   }
0101
0102   switch (prePoint->GetStepStatus())
0103   {
0104   case fPostStepDoItProc:
0105     if (savepoststepstatus != fGeomBoundary)
0106     {
0107       break;
0108     }
0109     else
0110     {
0111       std::cout << GetName() << ": New Hit for  " << std::endl;
0112       std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0113                 << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0114                 << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(saveprestepstatus)
0115                 << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(savepoststepstatus) << std::endl;
0116       std::cout << "last track: " << savetrackid
0117                 << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0118       std::cout << "phys pre vol: " << volume->GetName()
0119                 << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0120       std::cout << " previous phys pre vol: " << savevolpre->GetName()
0121                 << " previous phys post vol: " << savevolpost->GetName() << std::endl;
0122     }
0123     [[fallthrough]];
0124   case fGeomBoundary:
0125   case fUndefined:
0126     if (!hit)
0127     {
0128       hit = new PHG4Hitv1();
0129     }
0130     hit->set_layer(layer_id);
0131     // here we set the entrance values in cm
0132     hit->set_x(0, prePoint->GetPosition().x() / cm);
0133     hit->set_y(0, prePoint->GetPosition().y() / cm);
0134     hit->set_z(0, prePoint->GetPosition().z() / cm);
0135     // time in ns
0136     hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0137     // set the track ID
0138     hit->set_trkid(aTrack->GetTrackID());
0139     savetrackid = aTrack->GetTrackID();
0140     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0141     {
0142       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0143       {
0144         hit->set_trkid(pp->GetUserTrackId());
0145       }
0146     }
0147     // set the initial energy deposit
0148     edepsum = 0;
0149     if (whichactive > 0)
0150     {
0151       eionsum = 0;
0152       hit->set_eion(0);
0153       savehitcontainer = hits_;
0154     }
0155     else
0156     {
0157       std::cout << "implement stuff for whichactive < 0" << std::endl;
0158       gSystem->Exit(1);
0159     }
0160     if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0161     {
0162       if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0163       {
0164         hit->set_trkid(pp->GetUserTrackId());
0165         pp->GetShower()->add_g4hit_id(savehitcontainer->GetID(), hit->get_hit_id());
0166       }
0167     }
0168
0169     break;
0170   default:
0171     break;
0172   }
0173
0174   // some sanity checks for inconsistencies
0175   // check if this hit was created, if not print out last post step status
0176   if (!hit || !std::isfinite(hit->get_x(0)))
0177   {
0178     std::cout << GetName() << ": hit was not created" << std::endl;
0179     std::cout << "prestep status: " << PHG4StepStatusDecode::GetStepStatus(prePoint->GetStepStatus())
0180               << ", poststep status: " << PHG4StepStatusDecode::GetStepStatus(postPoint->GetStepStatus())
0181               << ", last pre step status: " << PHG4StepStatusDecode::GetStepStatus(saveprestepstatus)
0182               << ", last post step status: " << PHG4StepStatusDecode::GetStepStatus(savepoststepstatus) << std::endl;
0183     std::cout << "last track: " << savetrackid
0184               << ", current trackid: " << aTrack->GetTrackID() << std::endl;
0185     std::cout << "phys pre vol: " << volume->GetName()
0186               << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
0187     std::cout << " previous phys pre vol: " << savevolpre->GetName()
0188               << " previous phys post vol: " << savevolpost->GetName() << std::endl;
0189     gSystem->Exit(1);
0190   }
0191   // check if track id matches the initial one when the hit was created
0192   if (aTrack->GetTrackID() != savetrackid)
0193   {
0194     std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
0195     std::cout << "saved track: " << savetrackid
0196               << ", current trackid: " << aTrack->GetTrackID()
0197               << ", prestep status: " << prePoint->GetStepStatus()
0198               << ", previous post step status: " << savepoststepstatus
0199               << std::endl;
0200
0201     gSystem->Exit(1);
0202   }
0203   saveprestepstatus = prePoint->GetStepStatus();
0204   savepoststepstatus = postPoint->GetStepStatus();
0205   savevolpre = volume;
0206   savevolpost = touchpost->GetVolume();
0207
0208   // here we just update the exit values, it will be overwritten
0209   // for every step until we leave the volume or the particle
0210   // ceases to exist
0211   // sum up the energy to get total deposited
0212   edepsum += edep;
0213   if (whichactive > 0)
0214   {
0215     eionsum += eion;
0216   }
0217   // if any of these conditions is true this is the last step in
0218   // this volume and we need to save the hit
0219   // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
0220   // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
0221   // (happens when your detector goes outside world volume)
0222   // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
0223   // aTrack->GetTrackStatus() == fStopAndKill is also set)
0224   // aTrack->GetTrackStatus() == fStopAndKill: track ends
0225   if (postPoint->GetStepStatus() == fGeomBoundary ||
0226       postPoint->GetStepStatus() == fWorldBoundary ||
0227       postPoint->GetStepStatus() == fAtRestDoItProc ||
0228       aTrack->GetTrackStatus() == fStopAndKill)
0229   {
0230     // save only hits with energy deposit (or geantino)
0231     if (edepsum > 0 || geantino)
0232     {
0233       // update values at exit coordinates and set keep flag
0234       // of track to keep
0235       hit->set_x(1, postPoint->GetPosition().x() / cm);
0236       hit->set_y(1, postPoint->GetPosition().y() / cm);
0237       hit->set_z(1, postPoint->GetPosition().z() / cm);
0238
0239       hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0240       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0241       {
0242         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0243         {
0244           pp->SetKeep(1);  // we want to keep the track
0245         }
0246       }
0247       if (geantino)
0248       {
0249         hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0250         if (whichactive > 0)
0251         {
0252           hit->set_eion(-1);
0253         }
0254       }
0255       else
0256       {
0257         hit->set_edep(edepsum);
0258       }
0259       if (whichactive > 0)
0260       {
0261         hit->set_eion(eionsum);
0262       }
0263       savehitcontainer->AddHit(layer_id, hit);
0264       // ownership has been transferred to container, set to null
0265       // so we will create a new hit for the next track
0266       hit = nullptr;
0267     }
0268     else
0269     {
0270       // if this hit has no energy deposit, just reset it for reuse
0271       // this means we have to delete it in the dtor. If this was
0272       // the last hit we processed the memory is still allocated
0273       hit->Reset();
0274     }
0275   }
0276   // return true to indicate the hit was used
0277   return true;
0278 }
0279
0280 //____________________________________________________________________________..
0281 void PHG4PSTOFSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0282 {
0283   std::string hitnodename;
0284   if (detector_->SuperDetector() != "NONE")
0285   {
0286     hitnodename = "G4HIT_" + detector_->SuperDetector();
0287   }
0288   else
0289   {
0290     hitnodename = "G4HIT_" + detector_->GetName();
0291   }
0292
0293   // now look for the map and grab a pointer to it.
0294   hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0295
0296   // if we do not find the node we need to make it.
0297   if (!hits_)
0298   {
0299     std::cout << "PHG4PSTOFSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0300   }
0301 }