g4simulation/g4detectors/PHG4EnvelopeSteppingAction.cc

0001 #include "PHG4EnvelopeSteppingAction.h"
0002 #include "PHG4EnvelopeDetector.h"
0003
0004 #include <g4main/PHG4Hit.h>
0005 #include <g4main/PHG4HitContainer.h>
0006 #include <g4main/PHG4Hitv1.h>
0007 #include <g4main/PHG4Shower.h>
0008 #include <g4main/PHG4SteppingAction.h>  // for PHG4SteppingAction
0009
0010 #include <g4main/PHG4TrackUserInfoV1.h>
0011
0012 #include <phool/getClass.h>
0013
0014 #include <Geant4/G4ParticleDefinition.hh>      // for G4ParticleDefinition
0015 #include <Geant4/G4ReferenceCountedHandle.hh>  // for G4ReferenceCountedHandle
0016 #include <Geant4/G4Step.hh>
0017 #include <Geant4/G4StepPoint.hh>              // for G4StepPoint
0018 #include <Geant4/G4StepStatus.hh>             // for fGeomBoundary, fUndefined
0019 #include <Geant4/G4String.hh>                 // for G4String
0020 #include <Geant4/G4SystemOfUnits.hh>          // for mm, m
0021 #include <Geant4/G4ThreeVector.hh>            // for G4ThreeVector
0022 #include <Geant4/G4TouchableHandle.hh>        // for G4TouchableHandle
0023 #include <Geant4/G4Track.hh>                  // for G4Track
0024 #include <Geant4/G4TrackStatus.hh>            // for fStopAndKill
0025 #include <Geant4/G4Types.hh>                  // for G4double
0026 #include <Geant4/G4VTouchable.hh>             // for G4VTouchable
0027 #include <Geant4/G4VUserTrackInformation.hh>  // for G4VUserTrackInformation
0028
0029 #include <iostream>
0030 #include <string>  // for string, operator+, ope...
0031
0032 class G4VPhysicalVolume;
0033 class PHCompositeNode;
0034
0035 //______________________________________________________________
0036 PHG4EnvelopeSteppingAction::PHG4EnvelopeSteppingAction(PHG4EnvelopeDetector* detector)
0037   : PHG4SteppingAction(detector->GetName())
0038   , detector_(detector)
0039   , hits_(nullptr)
0040   , hit(nullptr)
0041 {
0042 }
0043
0044 //____________________________________________________________________________..
0045 bool PHG4EnvelopeSteppingAction::UserSteppingAction(const G4Step* aStep, bool /*was_used*/)
0046 {
0047   G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
0048   G4VPhysicalVolume* volume = touch->GetVolume();
0049
0050   int whichactive = detector_->IsInEnvelope(volume);
0051
0052   if (!whichactive)
0053   {
0054     return false;
0055   }
0056
0057   int layer_id = detector_->get_Layer();
0058   int tower_id = touch->GetCopyNumber();
0059
0060   /* Get energy deposited by this step */
0061   // G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
0062   G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
0063
0064   /* Get pointer to associated Geant4 track */
0065   const G4Track* aTrack = aStep->GetTrack();
0066
0067   // This detector is a black hole! Just put all kinetic energy into edep
0068   G4double edep = aTrack->GetKineticEnergy() / GeV;
0069   G4Track* killtrack = const_cast<G4Track*>(aTrack);
0070   killtrack->SetTrackStatus(fStopAndKill);
0071
0072   /* Make sure we are in a volume */
0073   if (detector_->IsActive())
0074   {
0075     /* Check if particle is 'geantino' */
0076     bool geantino = false;
0077     if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
0078     {
0079       geantino = true;
0080     }
0081
0082     /* Get Geant4 pre- and post-step points */
0083     G4StepPoint* prePoint = aStep->GetPreStepPoint();
0084     G4StepPoint* postPoint = aStep->GetPostStepPoint();
0085
0086     switch (prePoint->GetStepStatus())
0087     {
0088     case fGeomBoundary:
0089     case fUndefined:
0090       hit = new PHG4Hitv1();
0091       //      hit->set_layer(0);
0092       hit->set_scint_id(tower_id);
0093
0094       /* Set hit location (tower index) */
0095       // hit->set_index_j(idx_j);
0096       // hit->set_index_k(idx_k);
0097       // hit->set_index_l(idx_l);
0098       hit->set_index_j(0);
0099       hit->set_index_k(0);
0100       hit->set_index_l(0);
0101
0102       /* Set hit location (space point) */
0103       hit->set_x(0, prePoint->GetPosition().x() / cm);
0104       hit->set_y(0, prePoint->GetPosition().y() / cm);
0105       hit->set_z(0, prePoint->GetPosition().z() / cm);
0106
0107       /* Set momentum */
0108       hit->set_x(0, prePoint->GetMomentum().x() / GeV);
0109       hit->set_y(0, prePoint->GetMomentum().y() / GeV);
0110       hit->set_z(0, prePoint->GetMomentum().z() / GeV);
0111
0112       /* Set hit time */
0113       hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
0114
0115       /* set the track ID */
0116       {
0117         hit->set_trkid(aTrack->GetTrackID());
0118         if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0119         {
0120           if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0121           {
0122             hit->set_trkid(pp->GetUserTrackId());
0123             hit->set_shower_id(pp->GetShower()->get_id());
0124           }
0125         }
0126       }
0127
0128       /* set intial energy deposit */
0129       hit->set_edep(0);
0130       hit->set_eion(0);
0131
0132       hits_->AddHit(layer_id, hit);
0133
0134       {
0135         if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0136         {
0137           if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0138           {
0139             pp->GetShower()->add_g4hit_id(hits_->GetID(), hit->get_hit_id());
0140           }
0141         }
0142       }
0143
0144       break;
0145     default:
0146       break;
0147     }
0148
0149     /* Update exit values- will be overwritten with every step until
0150      * we leave the volume or the particle ceases to exist */
0151     hit->set_x(1, postPoint->GetPosition().x() / cm);
0152     hit->set_y(1, postPoint->GetPosition().y() / cm);
0153     hit->set_z(1, postPoint->GetPosition().z() / cm);
0154
0155     hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
0156
0157     /* sum up the energy to get total deposited */
0158     hit->set_edep(hit->get_edep() + edep);
0159     hit->set_eion(hit->get_eion() + eion);
0160
0161     if (geantino)
0162     {
0163       hit->set_edep(-1);  // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
0164       hit->set_eion(-1);
0165     }
0166     if (edep > 0)
0167     {
0168       if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
0169       {
0170         if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
0171         {
0172           pp->SetKeep(1);  // we want to keep the track
0173         }
0174       }
0175     }
0176     return true;
0177   }
0178   else
0179   {
0180     return false;
0181   }
0182 }
0183
0184 void PHG4EnvelopeSteppingAction::SetInterfacePointers(PHCompositeNode* topNode)
0185 {
0186   std::string hitnodename;
0187
0188   if (detector_->SuperDetector() != "NONE")
0189   {
0190     hitnodename = "G4HIT_ENVELOPE_" + detector_->SuperDetector();
0191   }
0192   else
0193   {
0194     hitnodename = "G4HIT_ENVELOPE_" + detector_->GetName();
0195   }
0196
0197   hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
0198
0199   // if we do not find the node it's messed up.
0200   if (!hits_)
0201   {
0202     std::cout << "PHG4CrystalCalorimeterSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
0203   }
0204 }