Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4EnvelopeSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4EnvelopeSteppingAction.cc
2 #include "PHG4EnvelopeDetector.h"
3 
4 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4Hitv1.h>
7 #include <g4main/PHG4Shower.h>
8 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
9 
11 
12 #include <phool/getClass.h>
13 
14 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
15 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
16 #include <Geant4/G4Step.hh>
17 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
18 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fUndefined
19 #include <Geant4/G4String.hh> // for G4String
20 #include <Geant4/G4SystemOfUnits.hh> // for mm, m
21 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
22 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
23 #include <Geant4/G4Track.hh> // for G4Track
24 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
25 #include <Geant4/G4Types.hh> // for G4double
26 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
27 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
28 
29 #include <iostream>
30 #include <string> // for string, operator+, ope...
31 
32 class G4VPhysicalVolume;
33 class PHCompositeNode;
34 
35 //______________________________________________________________
37  : PHG4SteppingAction(detector->GetName())
38  , detector_(detector)
39  , hits_(nullptr)
40  , hit(nullptr)
41 {
42 }
43 
44 //____________________________________________________________________________..
45 bool PHG4EnvelopeSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
46 {
47  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
48  G4VPhysicalVolume* volume = touch->GetVolume();
49 
50  int whichactive = detector_->IsInEnvelope(volume);
51 
52  if (!whichactive)
53  {
54  return false;
55  }
56 
57  int layer_id = detector_->get_Layer();
58  int tower_id = touch->GetCopyNumber();
59 
60  /* Get energy deposited by this step */
61  // G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
62  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
63 
64  /* Get pointer to associated Geant4 track */
65  const G4Track* aTrack = aStep->GetTrack();
66 
67  // This detector is a black hole! Just put all kinetic energy into edep
68  G4double edep = aTrack->GetKineticEnergy() / GeV;
69  G4Track* killtrack = const_cast<G4Track*>(aTrack);
70  killtrack->SetTrackStatus(fStopAndKill);
71 
72  /* Make sure we are in a volume */
73  if (detector_->IsActive())
74  {
75  /* Check if particle is 'geantino' */
76  bool geantino = false;
77  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
78  {
79  geantino = true;
80  }
81 
82  /* Get Geant4 pre- and post-step points */
83  G4StepPoint* prePoint = aStep->GetPreStepPoint();
84  G4StepPoint* postPoint = aStep->GetPostStepPoint();
85 
86  switch (prePoint->GetStepStatus())
87  {
88  case fGeomBoundary:
89  case fUndefined:
90  hit = new PHG4Hitv1();
91  // hit->set_layer(0);
92  hit->set_scint_id(tower_id);
93 
94  /* Set hit location (tower index) */
95  // hit->set_index_j(idx_j);
96  // hit->set_index_k(idx_k);
97  // hit->set_index_l(idx_l);
98  hit->set_index_j(0);
99  hit->set_index_k(0);
100  hit->set_index_l(0);
101 
102  /* Set hit location (space point) */
103  hit->set_x(0, prePoint->GetPosition().x() / cm);
104  hit->set_y(0, prePoint->GetPosition().y() / cm);
105  hit->set_z(0, prePoint->GetPosition().z() / cm);
106 
107  /* Set momentum */
108  hit->set_x(0, prePoint->GetMomentum().x() / GeV);
109  hit->set_y(0, prePoint->GetMomentum().y() / GeV);
110  hit->set_z(0, prePoint->GetMomentum().z() / GeV);
111 
112  /* Set hit time */
113  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
114 
115  /* set the track ID */
116  {
117  hit->set_trkid(aTrack->GetTrackID());
118  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
119  {
120  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
121  {
122  hit->set_trkid(pp->GetUserTrackId());
123  hit->set_shower_id(pp->GetShower()->get_id());
124  }
125  }
126  }
127 
128  /* set intial energy deposit */
129  hit->set_edep(0);
130  hit->set_eion(0);
131 
132  hits_->AddHit(layer_id, hit);
133 
134  {
135  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
136  {
137  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
138  {
139  pp->GetShower()->add_g4hit_id(hits_->GetID(), hit->get_hit_id());
140  }
141  }
142  }
143 
144  break;
145  default:
146  break;
147  }
148 
149  /* Update exit values- will be overwritten with every step until
150  * we leave the volume or the particle ceases to exist */
151  hit->set_x(1, postPoint->GetPosition().x() / cm);
152  hit->set_y(1, postPoint->GetPosition().y() / cm);
153  hit->set_z(1, postPoint->GetPosition().z() / cm);
154 
155  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
156 
157  /* sum up the energy to get total deposited */
158  hit->set_edep(hit->get_edep() + edep);
159  hit->set_eion(hit->get_eion() + eion);
160 
161  if (geantino)
162  {
163  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
164  hit->set_eion(-1);
165  }
166  if (edep > 0)
167  {
168  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
169  {
170  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
171  {
172  pp->SetKeep(1); // we want to keep the track
173  }
174  }
175  }
176  return true;
177  }
178  else
179  {
180  return false;
181  }
182 }
183 
185 {
186  std::string hitnodename;
187 
188  if (detector_->SuperDetector() != "NONE")
189  {
190  hitnodename = "G4HIT_ENVELOPE_" + detector_->SuperDetector();
191  }
192  else
193  {
194  hitnodename = "G4HIT_ENVELOPE_" + detector_->GetName();
195  }
196 
197  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
198 
199  // if we do not find the node it's messed up.
200  if (!hits_)
201  {
202  std::cout << "PHG4CrystalCalorimeterSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
203  }
204 }