Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4SectorSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4SectorSteppingAction.cc
2 #include "PHG4SectorDetector.h"
3 
4 #include <g4main/PHG4Hit.h>
6 #include <g4main/PHG4Hitv1.h>
7 #include <g4main/PHG4Shower.h>
8 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
10 
11 #include <phool/getClass.h>
12 
13 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
14 #include <Geant4/G4Step.hh>
15 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
16 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRestD...
17 #include <Geant4/G4String.hh> // for G4String
18 #include <Geant4/G4SystemOfUnits.hh> // for cm, GeV, nanosecond
19 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
20 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
21 #include <Geant4/G4Track.hh> // for G4Track
22 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
23 #include <Geant4/G4Types.hh> // for G4double
24 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
25 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
26 
27 #include <iostream>
28 #include <string> // for string, operator+, oper...
29 
30 class G4VPhysicalVolume;
31 class PHCompositeNode;
32 
33 //____________________________________________________________________________..
35  : PHG4SteppingAction(detector->GetName())
36  , detector_(detector)
37 {
38 }
39 
41 {
42  // if the last hit was a zero energie deposit hit, it is just reset
43  // and the memory is still allocated, so we need to delete it here
44  // if the last hit was saved, hit is a nullptr pointer which are
45  // legal to delete (it results in a no operation)
46  delete hit;
47 }
48 
49 //____________________________________________________________________________..
50 bool PHG4SectorSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
51 {
52  // get volume of the current step
53  G4VPhysicalVolume* volume =
54  aStep->GetPreStepPoint()->GetTouchableHandle()->GetVolume();
55 
56  // collect energy and track length step by step
57  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
58  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
59 
60  const G4Track* aTrack = aStep->GetTrack();
61 
62  // make sure we are in a volume
63  if (detector_->IsInSectorActive(volume))
64  {
65  bool geantino = false;
66  // the check for the pdg code speeds things up, I do not want to make
67  // an expensive string compare for every track when we know
68  // geantino or chargedgeantino has pid=0
69  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
70  {
71  geantino = true;
72  }
73  G4StepPoint* prePoint = aStep->GetPreStepPoint();
74  G4StepPoint* postPoint = aStep->GetPostStepPoint();
75  // std::cout << "track id " << aTrack->GetTrackID() << std::endl;
76  // std::cout << "time prepoint: " << prePoint->GetGlobalTime() << std::endl;
77  // std::cout << "time postpoint: " << postPoint->GetGlobalTime() << std::endl;
78  //layer_id is sector number
79  switch (prePoint->GetStepStatus())
80  {
81  case fGeomBoundary:
82  case fUndefined:
83  if (!hit)
84  {
85  hit = new PHG4Hitv1();
86  }
87  //here we set the entrance values in cm
88  hit->set_x(0, prePoint->GetPosition().x() / cm);
89  hit->set_y(0, prePoint->GetPosition().y() / cm);
90  hit->set_z(0, prePoint->GetPosition().z() / cm);
91  // time in ns
92  hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
93  //set the track ID
94  hit->set_trkid(aTrack->GetTrackID());
95  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
96  {
97  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
98  {
99  hit->set_trkid(pp->GetUserTrackId());
100  hit->set_shower_id(pp->GetShower()->get_id());
101  saveshower = pp->GetShower();
102  }
103  }
104 
105  //set the initial energy deposit
106  hit->set_edep(0);
107  hit->set_eion(0); // only implemented for v5 otherwise empty
108  layer_id = aStep->GetPreStepPoint()->GetTouchable()->GetReplicaNumber(1);
109  // hit->set_light_yield(0);
110 
111  break;
112  default:
113  break;
114  }
115  // here we just update the exit values, it will be overwritten
116  // for every step until we leave the volume or the particle
117  // ceases to exist
118  hit->set_x(1, postPoint->GetPosition().x() / cm);
119  hit->set_y(1, postPoint->GetPosition().y() / cm);
120  hit->set_z(1, postPoint->GetPosition().z() / cm);
121 
122  hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
123  //sum up the energy to get total deposited
124  hit->set_edep(hit->get_edep() + edep);
125  hit->set_eion(hit->get_eion() + eion);
126  hit->set_path_length(aTrack->GetTrackLength() / cm);
127  if (geantino)
128  {
129  hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
130  }
131  if (edep > 0)
132  {
133  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
134  {
135  if (PHG4TrackUserInfoV1* pp =
136  dynamic_cast<PHG4TrackUserInfoV1*>(p))
137  {
138  pp->SetKeep(1); // we want to keep the track
139  }
140  }
141  }
142  // if any of these conditions is true this is the last step in
143  // this volume and we need to save the hit
144  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
145  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
146  // (not sure if this will ever be the case)
147  // aTrack->GetTrackStatus() == fStopAndKill: track ends
148  if (postPoint->GetStepStatus() == fGeomBoundary ||
149  postPoint->GetStepStatus() == fWorldBoundary ||
150  postPoint->GetStepStatus() == fAtRestDoItProc ||
151  aTrack->GetTrackStatus() == fStopAndKill)
152  {
153  // save only hits with energy deposit (or -1 for geantino)
154  if (hit->get_edep())
155  {
157  if (saveshower)
158  {
160  }
161  // ownership has been transferred to container, set to null
162  // so we will create a new hit for the next track
163  hit = nullptr;
164  }
165  else
166  {
167  // if this hit has no energy deposit, just reset it for reuse
168  // this means we have to delete it in the dtor. If this was
169  // the last hit we processed the memory is still allocated
170  hit->Reset();
171  }
172  }
173 
174  // hit->identify();
175  // return true to indicate the hit was used
176  return true;
177  }
178  else
179  {
180  return false;
181  }
182 }
183 
184 //____________________________________________________________________________..
186 {
187  std::string hitnodename;
188  if (detector_->SuperDetector() != "NONE")
189  {
190  hitnodename = "G4HIT_" + detector_->SuperDetector();
191  }
192  else
193  {
194  hitnodename = "G4HIT_" + detector_->GetName();
195  }
196 
197  //now look for the map and grab a pointer to it.
198  hits_ = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
199 
200  // if we do not find the node we need to make it.
201  if (!hits_)
202  {
203  std::cout << "PHG4SectorSteppingAction::SetTopNode - unable to find "
204  << hitnodename << std::endl;
205  }
206 }