ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4InttSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4InttSteppingAction.cc
2 #include "PHG4InttDefs.h"
3 #include "PHG4InttDetector.h"
4 
5 #include <phparameter/PHParameters.h>
6 #include <phparameter/PHParametersContainer.h>
7 
8 #include <g4main/PHG4Hit.h>
10 #include <g4main/PHG4Hitv1.h>
11 #include <g4main/PHG4Shower.h>
12 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
14 
15 #include <phool/getClass.h>
16 
17 #include <TSystem.h>
18 
19 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
20 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
21 #include <Geant4/G4Step.hh>
22 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
23 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRes...
24 #include <Geant4/G4String.hh> // for G4String
25 #include <Geant4/G4SystemOfUnits.hh>
26 #include <Geant4/G4ThreeVector.hh>
27 #include <Geant4/G4TouchableHandle.hh>
28 #include <Geant4/G4Track.hh> // for G4Track
29 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
30 #include <Geant4/G4Types.hh> // for G4double
31 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
32 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
33 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
34 
35 #include <cassert> // for assert
36 #include <iostream>
37 #include <set> // for set
38 #include <string> // for operator<<, string
39 #include <tuple> // for tie, tuple
40 
41 class PHCompositeNode;
42 
43 using namespace std;
44 
45 //____________________________________________________________________________..
46 PHG4InttSteppingAction::PHG4InttSteppingAction(PHG4InttDetector* detector, const PHParametersContainer* parameters, const pair<vector<pair<int, int>>::const_iterator, vector<pair<int, int>>::const_iterator>& layer_begin_end)
47  : PHG4SteppingAction(detector->GetName())
48  , m_Detector(detector)
49  , m_Hits(nullptr)
50  , m_AbsorberHits(nullptr)
51  , m_Hit(nullptr)
52  , m_SaveHitContainer(nullptr)
53  , m_SaveShower(nullptr)
54  , m_ParamsContainer(parameters)
55 {
56  // loop over layers to get laddertype nd active status for each layer
57  for (auto layeriter = layer_begin_end.first; layeriter != layer_begin_end.second; ++layeriter)
58  {
59  int layer = layeriter->second;
60  const PHParameters* par = m_ParamsContainer->GetParameters(layer);
61  m_IsActiveMap[layer] = par->get_int_param("active");
62  m_IsBlackHoleMap[layer] = par->get_int_param("blackhole");
63  m_LadderTypeMap.insert(make_pair(layer, par->get_int_param("laddertype")));
64  m_InttToTrackerLayerMap.insert(make_pair(layeriter->second, layeriter->first));
65  }
66 
67  // Get the parameters for each laddertype
68  for (auto iter = PHG4InttDefs::m_SensorSegmentationSet.begin(); iter != PHG4InttDefs::m_SensorSegmentationSet.end(); ++iter)
69  {
70  const PHParameters* par = m_ParamsContainer->GetParameters(*iter);
71  m_StripYMap.insert(make_pair(*iter, par->get_double_param("strip_y") * cm));
72  m_StripZMap.insert(make_pair(*iter, make_pair(par->get_double_param("strip_z_0") * cm, par->get_double_param("strip_z_1") * cm)));
73  m_nStripsPhiCell.insert(make_pair(*iter, par->get_int_param("nstrips_phi_cell")));
74  m_nStripsZSensor.insert(make_pair(*iter, make_pair(par->get_int_param("nstrips_z_sensor_0"), par->get_int_param("nstrips_z_sensor_1"))));
75  }
76 }
77 
79 {
80  // if the last hit was a zero energie deposit hit, it is just reset
81  // and the memory is still allocated, so we need to delete it here
82  // if the last hit was saved, hit is a nullptr pointer which are
83  // legal to delete (it results in a no operation)
84  delete m_Hit;
85 }
86 
87 //____________________________________________________________________________..
89 {
91  // get volume of the current step
93  const G4Track* aTrack = aStep->GetTrack();
94  G4StepPoint* prePoint = aStep->GetPreStepPoint();
95  G4StepPoint* postPoint = aStep->GetPostStepPoint();
96 
97  const int whichactive = m_Detector->IsInIntt(volume);
98 
99  if (!whichactive)
100  {
101  return false;
102  }
103 
104  if (Verbosity() > 0)
105  {
106  cout << endl
107  << "PHG4SilicoTrackerSteppingAction::UserSteppingAction for volume name (pre) " << touch->GetVolume()->GetName()
108  << " volume name (1) " << touch->GetVolume(1)->GetName()
109  << " volume->GetTranslation " << touch->GetVolume()->GetTranslation()
110  << endl;
111  }
112 
113  // set ladder index
114  int sphxlayer = 0;
115  int inttlayer = 0;
116  int ladderz = 0;
117  int ladderphi = 0;
118  int zposneg = 0;
119 
120  if (whichactive > 0) // silicon active sensor
121  {
122  // Get the layer and ladder information which is step up in the volume hierarchy
123  // the ladder also contains inactive volumes but we check in m_Detector->IsInIntt(volume)
124  // if we are in an active logical volume whioch is located in this ladder
125  auto iter = m_Detector->get_ActiveVolumeTuple(touch->GetVolume(1));
126  tie(inttlayer, ladderz, ladderphi, zposneg) = iter->second;
127  if (Verbosity() > 0)
128  cout << " inttlayer " << inttlayer << " ladderz_base " << ladderz << " ladderphi " << ladderphi << " zposneg " << zposneg << endl;
129  if (inttlayer < 0 || inttlayer > 7)
130  {
131  assert(!"PHG4InttSteppingAction: check Intt ladder layer.");
132  }
133  sphxlayer = m_InttToTrackerLayerMap.find(inttlayer)->second;
134  map<int, int>::const_iterator activeiter = m_IsActiveMap.find(inttlayer);
135  if (activeiter == m_IsActiveMap.end())
136  {
137  cout << "PHG4InttSteppingAction: could not find active flag for layer " << inttlayer << endl;
138  gSystem->Exit(1);
139  }
140  if (activeiter->second == 0)
141  {
142  return false;
143  }
144  }
145  else // whichactive < 0, silicon inactive area, FPHX, stabe etc. as absorbers
146  {
147  auto iter = m_Detector->get_PassiveVolumeTuple(touch->GetVolume(0)->GetLogicalVolume());
148  tie(inttlayer, ladderz) = iter->second;
149  sphxlayer = inttlayer; //for absorber we use the Intt layer, not the tracking layer in sPHENIX
150  } // end of si inactive area block
151 
152  // collect energy and track length step by step
153  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
154  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
155 
156  // if this block stops everything, just put all kinetic energy into edep
157  if ((m_IsBlackHoleMap.find(inttlayer))->second == 1)
158  {
159  edep = aTrack->GetKineticEnergy() / GeV;
160  G4Track* killtrack = const_cast<G4Track*>(aTrack);
161  killtrack->SetTrackStatus(fStopAndKill);
162  }
163 
164  bool geantino = false;
165 
166  // the check for the pdg code speeds things up, I do not want to make
167  // an expensive string compare for every track when we know
168  // geantino or chargedgeantino has pid=0
169  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 && aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
170  {
171  geantino = true;
172  }
173 
174  if (Verbosity() > 1)
175  {
176  cout << " aTrack->GetTrackID " << aTrack->GetTrackID() << " aTrack->GetParentID " << aTrack->GetParentID()
177  << " Intt layer " << inttlayer
178  << " prePoint step status = " << prePoint->GetStepStatus() << " postPoint step status = " << postPoint->GetStepStatus()
179  << endl;
180  }
181  switch (prePoint->GetStepStatus())
182  {
183  case fGeomBoundary:
184  case fUndefined:
185 
186  if (Verbosity() > 1)
187  {
188  cout << " start new g4hit for: aTrack->GetParentID " << aTrack->GetParentID() << " aTrack->GetTrackID " << aTrack->GetTrackID() << " Intt layer " << inttlayer
189  << " prePoint step status = " << prePoint->GetStepStatus() << " postPoint->GetStepStatus = " << postPoint->GetStepStatus() << endl;
190  }
191  // if previous hit was saved, hit pointer was set to nullptr
192  // and we have to make a new one
193  if (!m_Hit)
194  {
195  m_Hit = new PHG4Hitv1();
196  }
197 
198  // set the index values needed to locate the sensor strip
199  if (zposneg == 1)
200  {
201  ladderz += 2; // ladderz = 0, 1 for negative z and = 2, 3 for positive z
202  }
203  if(Verbosity() > 0) cout << " ladderz = " << ladderz << endl;
204 
205  m_Hit->set_ladder_z_index(ladderz);
206 
207  if (whichactive > 0)
208  {
209  m_Hit->set_layer(sphxlayer);
210  m_Hit->set_ladder_phi_index(ladderphi);
211  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
212  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
213  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
214  m_Hit->set_eion(0);
215  }
216 
217  //here we set the entrance values in cm
218  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
219  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
220  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
221 
222  StoreLocalCoordinate(m_Hit, aStep, true, false);
223 
224  if (Verbosity() > 0)
225  {
226  cout << " prePoint hit position x,y,z = " << prePoint->GetPosition().x() / cm
227  << " " << prePoint->GetPosition().y() / cm
228  << " " << prePoint->GetPosition().z() / cm
229  << endl;
230  }
231 
232  // time in ns
233  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
234 
235  //set the track ID
236  m_Hit->set_trkid(aTrack->GetTrackID());
237 
238  //set the initial energy deposit
239  m_Hit->set_edep(0);
240 
241  if (whichactive > 0) // return of IsInIntt, > 0 hit in si-strip, < 0 hit in absorber
242  {
243  // Now save the container we want to add this hit to
245  }
246  else
247  {
249  }
250 
252  {
253  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
254  {
255  m_Hit->set_trkid(pp->GetUserTrackId());
256  m_Hit->set_shower_id(pp->GetShower()->get_id());
257  m_SaveShower = pp->GetShower();
258  }
259  }
260  break;
261 
262  default:
263  break;
264  }
265 
266  //cout << " Update exit values for prePoint->GetStepStatus = " << prePoint->GetStepStatus() << " and postPoint->GetStepStatus = " << postPoint->GetStepStatus() << endl;
267 
268  // here we just update the exit values, it will be overwritten
269  // for every step until we leave the volume or the particle
270  // ceases to exist
271  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
272  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
273  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
274 
275  StoreLocalCoordinate(m_Hit, aStep, false, true);
276 
277  if (whichactive > 0)
278  {
279  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
280  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
281  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
282  m_Hit->set_eion(m_Hit->get_eion() + eion);
283  }
284 
285  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
286 
287  //sum up the energy to get total deposited
289 
290  if (geantino)
291  {
292  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
293  if (whichactive > 0)
294  {
295  m_Hit->set_eion(-1);
296  }
297  }
298 
299  if (edep > 0)
300  {
302  {
303  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
304  {
305  pp->SetKeep(1); // we want to keep the track
306  }
307  }
308  }
309 
310  // if any of these conditions is true this is the last step in
311  // this volume and we need to save the hit
312  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
313  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
314  // (happens when your detector goes outside world volume)
315  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
316  // aTrack->GetTrackStatus() == fStopAndKill is also set)
317  // aTrack->GetTrackStatus() == fStopAndKill: track ends
318  if (postPoint->GetStepStatus() == fGeomBoundary ||
319  postPoint->GetStepStatus() == fWorldBoundary ||
320  postPoint->GetStepStatus() == fAtRestDoItProc ||
321  aTrack->GetTrackStatus() == fStopAndKill)
322  {
323  if (Verbosity() > 0)
324  {
325  cout << " postPoint step status changed to " << postPoint->GetStepStatus() << " or aTrack status changed to "
326  << aTrack->GetTrackStatus() << endl;
327  cout << " end g4hit for: aTrack->GetParentID " << aTrack->GetParentID() << " aTrack->GetTrackID " << aTrack->GetTrackID() << " eion = " << eion << endl;
328  cout << " end hit position x,y,z = " << postPoint->GetPosition().x() / cm
329  << " " << postPoint->GetPosition().y() / cm
330  << " " << postPoint->GetPosition().z() / cm
331  << endl;
332  }
333 
334  // save only hits with energy deposit (or -1 for geantino)
335  if (m_Hit->get_edep())
336  {
337  m_SaveHitContainer->AddHit(sphxlayer, m_Hit);
338  if (m_SaveShower)
339  {
341  }
342  if (Verbosity() > 0)
343  {
344  m_Hit->print();
345  }
346  // ownership has been transferred to container, set to null
347  // so we will create a new hit for the next track
348  m_Hit = nullptr;
349  }
350  else
351  {
352  // if this hit has no energy deposit, just reset it for reuse
353  // this means we have to delete it in the dtor. If this was
354  // the last hit we processed the memory is still allocated
355  m_Hit->Reset();
356  }
357  }
358  return true;
359 }
360 
361 //____________________________________________________________________________..
363 {
364  const string detectorname = (m_Detector->SuperDetector() != "NONE") ? m_Detector->SuperDetector() : m_Detector->GetName();
365  const string hitnodename = "G4HIT_" + detectorname;
366  const string absorbernodename = "G4HIT_ABSORBER_" + detectorname;
367 
368  //now look for the map and grab a pointer to it.
369  m_Hits = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
370  m_AbsorberHits = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
371 
372  // if we do not find the node it's messed up.
373  if (!m_Hits)
374  cout << "PHG4InttSteppingAction::SetTopNode - unable to find " << hitnodename << endl;
375 
376  if (!m_AbsorberHits && Verbosity() > 1)
377  cout << "PHG4InttSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
378 }