ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4MvtxSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4MvtxSteppingAction.cc
2 
3 #include "PHG4MvtxDetector.h"
4 
5 #include <g4main/PHG4Hit.h>
7 #include <g4main/PHG4Hitv1.h>
8 #include <g4main/PHG4Shower.h>
9 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
11 
12 #include <phool/getClass.h>
13 #include <phool/phool.h> // for PHWHERE
14 
15 #include <Geant4/G4NavigationHistory.hh>
16 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
17 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
18 #include <Geant4/G4Step.hh>
19 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
20 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
21 #include <Geant4/G4String.hh> // for G4String
22 #include <Geant4/G4SystemOfUnits.hh>
23 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
24 #include <Geant4/G4TouchableHandle.hh> // for G4TouchableHandle
25 #include <Geant4/G4Track.hh> // for G4Track
26 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
27 #include <Geant4/G4Types.hh> // for G4double
28 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
29 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
30 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
31 
32 #include <boost/tokenizer.hpp>
33 // this is an ugly hack, the gcc optimizer has a bug which
34 // triggers the uninitialized variable warning which
35 // stops compilation because of our -Werror
36 #include <boost/version.hpp> // to get BOOST_VERSION
37 #if (__GNUC__ == 4 && __GNUC_MINOR__ == 4 && BOOST_VERSION == 105700)
38 #pragma GCC diagnostic ignored "-Wuninitialized"
39 #pragma message "ignoring bogus gcc warning in boost header lexical_cast.hpp"
40 #include <boost/lexical_cast.hpp>
41 #pragma GCC diagnostic warning "-Wuninitialized"
42 #else
43 #include <boost/lexical_cast.hpp>
44 #endif
45 
46 #include <cmath> // for NAN
47 #include <cstdlib> // for exit
48 #include <iostream>
49 #include <string> // for operator<<, basic_string
50 
51 class PHCompositeNode;
52 
53 using namespace std;
54 //____________________________________________________________________________..
56  : PHG4SteppingAction(detector->GetName())
57  , m_Detector(detector)
58  , m_HitContainer(nullptr)
59  , m_AbsorberhitContainer(nullptr)
60  , m_Hit(nullptr)
61  , m_SaveShower(nullptr)
62 {
63 }
64 
65 //____________________________________________________________________________..
67 {
68  // if the last hit was a zero energie deposit hit, it is just reset
69  // and the memory is still allocated, so we need to delete it here
70  // if the last hit was saved, hit is a nullptr pointer which are
71  // legal to delete (it results in a no operation)
72  delete m_Hit;
73 }
74 
75 //____________________________________________________________________________..
77 {
79  // get volume of the current step
80  G4VPhysicalVolume* sensor_volume = touch->GetVolume();
81 
82  // PHG4MvtxDetector->IsInMvtx(volume)
83  // returns
84  // 0 if outside of Mvtx
85  // 1 if inside sensor
86 
87  // This checks if the volume is a sensor (doesn't tell us unique layer)
88  // PHG4MvtxDetector->IsSensor(volume)
89  // returns
90  // 1 if volume is a sensor
91  // 0 if not
92  int whichactive = m_Detector->IsSensor(sensor_volume);
93 
94  if (!whichactive)
95  {
96  return false;
97  }
98 
99  // This tells us if the volume belongs to the right stave for this layer
100  // From the GDML file the 3rd volume up should be the half-stave
101  // PHG4MvtxTelescopeDetector_->IsInMvtx(volume)
102  // returns
103  // 1 if in ladder belonging to this layer
104  // 0 if not
105  int layer_id = -9999;
106  int stave_id = -9999;
107  //cout << endl << " In UserSteppingAction for layer " << layer_id << endl;
108  G4VPhysicalVolume* vstave = touch->GetVolume(3);
109  whichactive = m_Detector->IsInMvtx(vstave, layer_id, stave_id);
110  if (layer_id < 0 || stave_id < 0)
111  {
112  cout << PHWHERE << "invalid Mvtx's layer (" << layer_id << ") or stave (" << stave_id << ") index " << endl;
113  exit(1);
114  }
115 
116  if (!whichactive)
117  {
118  return false;
119  }
120 
121  if (Verbosity() > 5)
122  {
123  // make sure we know where we are!
124  G4VPhysicalVolume* vtest = touch->GetVolume();
125  cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume " << vtest->GetName() << endl;
126  G4VPhysicalVolume* vtest1 = touch->GetVolume(1);
127  cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 1 up " << vtest1->GetName() << endl;
128  G4VPhysicalVolume* vtest2 = touch->GetVolume(2);
129  cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 2 up " << vtest2->GetName() << endl;
130  G4VPhysicalVolume* vtest3 = touch->GetVolume(3);
131  cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 3 up " << vtest3->GetName() << endl;
132  G4VPhysicalVolume* vtest4 = touch->GetVolume(4);
133  cout << "Entering PHG4MvtxSteppingAction::UserSteppingAction for volume 4 up " << vtest4->GetName() << endl;
134  }
135 
136  //=======================================================================
137  // We want the location of the hit
138  // Here we will record in the hit object:
139  // The stave, half-stave, module and chip numbers
140  // The energy deposited
141  // The entry point and exit point in world coordinates
142  // The entry point and exit point in local (sensor) coordinates
143  // The pixel number will be derived later from the entry and exit points in the sensor local coordinates
144  //=======================================================================
145 
146  //int stave_number = -1; // stave_number from stave map values
147  int half_stave_number = -1;
148  int module_number = -1;
149  int chip_number = -1;
150 
151  if (Verbosity() > 0)
152  cout << endl
153  << " UserSteppingAction: layer " << layer_id;
154  boost::char_separator<char> sep("_");
155  boost::tokenizer<boost::char_separator<char> >::const_iterator tokeniter;
156 
157  //OLD ITS.gdml: chip number is from "ITSUChip[layer number]_[chip number]
158  //NEW: chip number is from "MVTXChip_[chip number]
159  G4VPhysicalVolume* v1 = touch->GetVolume(1);
160  boost::tokenizer<boost::char_separator<char> > tok1(v1->GetName(), sep);
161  tokeniter = tok1.begin();
162  ++tokeniter;
163  chip_number = boost::lexical_cast<int>(*tokeniter);
164  if (Verbosity() > 0) cout << " chip " << chip_number;
165  G4VPhysicalVolume* v2 = touch->GetVolume(2);
166  boost::tokenizer<boost::char_separator<char> > tok2(v2->GetName(), sep);
167  tokeniter = tok2.begin();
168  ++tokeniter;
169  module_number = boost::lexical_cast<int>(*tokeniter);
170  if (Verbosity() > 0) cout << " module " << module_number;
171 
172  // The stave number is the imprint number from the assembly volume imprint
173  // The assembly volume history string format is (e.g.):
174  // av_13_impr_4_ITSUHalfStave6_pv_1
175  // where "impr_4" says stave number 4
176  // and where "ITSUHalfStave6_pv_1" says hald stave number 1 in layer number 6
177 
178  G4VPhysicalVolume* v3 = touch->GetVolume(3);
179  boost::tokenizer<boost::char_separator<char> > tok3(v3->GetName(), sep);
180  tokeniter = tok3.begin();
181  ++tokeniter;
182  ++tokeniter;
183  ++tokeniter;
184  //stave_number = boost::lexical_cast<int>(*tokeniter) - 1; // starts counting imprints at 1, we count staves from 0!
185  if (Verbosity() > 0) cout << " stave " << stave_id;
186  ++tokeniter;
187  ++tokeniter;
188  ++tokeniter;
189  half_stave_number = boost::lexical_cast<int>(*tokeniter);
190  if (Verbosity() > 0) cout << " half_stave " << half_stave_number;
191 
192  // FYI: doing string compares inside a stepping action sounds like a recipe
193  // for failure inside a heavy ion event... we'll wait and see how badly
194  // this profiles. -MPM
195 
196  // Now we want to collect information about the hit
197 
198  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
199  const G4Track* aTrack = aStep->GetTrack();
200  if (Verbosity() > 0) cout << " edep = " << edep << endl;
201 
202  // if this cylinder stops everything, just put all kinetic energy into edep
203  if (m_Detector->IsBlackHole(layer_id))
204  {
205  edep = aTrack->GetKineticEnergy() / GeV;
206  G4Track* killtrack = const_cast<G4Track*>(aTrack);
207  killtrack->SetTrackStatus(fStopAndKill);
208  }
209 
210  // test if we are active
211  if (m_Detector->IsActive(layer_id))
212  {
213  bool geantino = false;
214  // the check for the pdg code speeds things up, I do not want to make
215  // an expensive string compare for every track when we know
216  // geantino or chargedgeantino has pid=0
217  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
218  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != string::npos)
219  {
220  geantino = true;
221  }
222 
223  G4ThreeVector worldPosition; // localPosition;
224  G4TouchableHandle theTouchable;
225 
226  G4StepPoint* prePoint = aStep->GetPreStepPoint();
227  G4StepPoint* postPoint = aStep->GetPostStepPoint();
228 
229  if (Verbosity() > 0)
230  {
231  G4ParticleDefinition* def = aTrack->GetDefinition();
232  cout << " Particle: " << def->GetParticleName() << endl;
233  }
234 
235  switch (prePoint->GetStepStatus())
236  {
237  case fGeomBoundary:
238  case fUndefined:
239  if (!m_Hit)
240  {
241  m_Hit = new PHG4Hitv1();
242  }
243  m_Hit->set_layer((unsigned int) layer_id);
244 
245  // set the index values needed to locate the sensor strip
250 
251  worldPosition = prePoint->GetPosition();
252 
253  if (Verbosity() > 0)
254  {
255  theTouchable = prePoint->GetTouchableHandle();
256  cout << "entering: depth = " << theTouchable->GetHistory()->GetDepth() << endl;
257  G4VPhysicalVolume *vol1 = theTouchable->GetVolume();
258  cout << "entering volume name = " << vol1->GetName() << endl;
259  }
260 
261  /*
262  localPosition = theTouchable->GetHistory()->GetTopTransform().TransformPoint(worldPosition);
263  */
264 
265  // Store the local coordinates for the entry point
266  StoreLocalCoordinate(m_Hit, aStep, true, false);
267 
268  // Store the entrance values in cm in world coordinates
269  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
270  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
271  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
272 
273  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
274  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
275  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
276 
277  // time in ns
278  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
279  //set the track ID
280  m_Hit->set_trkid(aTrack->GetTrackID());
282  {
283  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
284  {
285  m_Hit->set_trkid(pp->GetUserTrackId());
286  m_Hit->set_shower_id(pp->GetShower()->get_id());
287  m_SaveShower = pp->GetShower();
288  }
289  }
290  //set the initial energy deposit
291  m_Hit->set_edep(0);
292 
293  // Now add the hit
294  // m_Hit->print();
295  //m_HitContainer->AddHit(layer_id, m_Hit);
296  break;
297  default:
298  break;
299  }
300 
301  // here we just update the exit values, it will be overwritten
302  // for every step until we leave the volume or the particle
303  // ceases to exist
304 
305  /*
306  // Note that the after you reach the boundary the touchable for the postPoint points to the next volume, not the sensor!
307  // This was given to me as the way to get back to the sensor volume, but it does not work
308  theTouchable = postPoint->GetTouchableHandle();
309  localPosition = history->GetTransform(history->GetDepth() - 1).TransformPoint(worldPosition);
310  cout << "Exit local coords: x " << localPosition.x() / cm << " y " << localPosition.y() / cm << " z " << localPosition.z() / cm << endl;
311  */
312 
313  /*
314  // Use the prePoint from the final step for now, until I understand how to get the exit point in the sensor volume coordinates
315  //======================================================================================
316  theTouchable = prePoint->GetTouchableHandle();
317  vol2 = theTouchable->GetVolume();
318  if(Verbosity() > 0)
319  cout << "exiting volume name = " << vol2->GetName() << endl;
320  worldPosition = prePoint->GetPosition();
321  if(Verbosity() > 0)
322  cout << "Exit world coords prePoint: x " << worldPosition.x() / cm << " y " << worldPosition.y() / cm << " z " << worldPosition.z() / cm << endl;
323 
324  // This is for consistency with the local coord position, the world coordinate exit position is correct
325  m_Hit->set_x( 1, prePoint->GetPosition().x() / cm );
326  m_Hit->set_y( 1, prePoint->GetPosition().y() / cm );
327  m_Hit->set_z( 1, prePoint->GetPosition().z() / cm );
328 
329  const G4NavigationHistory *history = theTouchable->GetHistory();
330  //cout << "exiting: depth = " << history->GetDepth() << " volume name = " << history->GetVolume(history->GetDepth())->GetName() << endl;
331  localPosition = history->GetTransform(history->GetDepth()).TransformPoint(worldPosition);
332 
333  m_Hit->set_local_x(1, localPosition.x() / cm);
334  m_Hit->set_local_y(1, localPosition.y() / cm);
335  m_Hit->set_local_z(1, localPosition.z() / cm);
336  */
337 
338  // Store the local coordinates for the exit point
339  StoreLocalCoordinate(m_Hit, aStep, false, true);
340 
341  // Store world coordinates for the exit point
342  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
343  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
344  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
345 
346  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
347  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
348  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
349 
350  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
351  //sum up the energy to get total deposited
353  if (geantino)
354  {
355  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
356  }
357  if (edep > 0)
358  {
360  {
361  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
362  {
363  pp->SetKeep(1); // we want to keep the track
364  }
365  }
366  }
367 
368  if (Verbosity() > 0)
369  {
370  G4StepPoint* prePoint = aStep->GetPreStepPoint();
371  G4StepPoint* postPoint = aStep->GetPostStepPoint();
372  cout << "----- PHg4MvtxSteppingAction::UserSteppingAction - active volume = " << sensor_volume->GetName() << endl;
373  cout << " layer = " << layer_id << endl;
374  cout << " stave number = " << stave_id << " half_stave_number = " << half_stave_number << endl;
375  cout << " module number = " << module_number << endl;
376  cout << " chip number = " << chip_number << endl;
377  cout << " prepoint x position " << prePoint->GetPosition().x() / cm << endl;
378  cout << " prepoint y position " << prePoint->GetPosition().y() / cm << endl;
379  cout << " prepoint z position " << prePoint->GetPosition().z() / cm << endl;
380  cout << " postpoint x position " << postPoint->GetPosition().x() / cm << endl;
381  cout << " postpoint y position " << postPoint->GetPosition().y() / cm << endl;
382  cout << " postpoint z position " << postPoint->GetPosition().z() / cm << endl;
383  cout << " edep " << edep << endl;
384  }
385 
386  if (Verbosity() > 0)
387  {
388  cout << " stepping action found hit:" << endl;
389  m_Hit->print();
390  cout << endl
391  << endl;
392  }
393 
394  // if any of these conditions is true this is the last step in
395  // this volume and we need to save the hit
396  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
397  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
398  // (happens when your detector goes outside world volume)
399  // postPoint->GetStepStatus() == fAtRestDoItProc: track stops (typically
400  // aTrack->GetTrackStatus() == fStopAndKill is also set)
401  // aTrack->GetTrackStatus() == fStopAndKill: track ends
402  if (postPoint->GetStepStatus() == fGeomBoundary ||
403  postPoint->GetStepStatus() == fWorldBoundary ||
404  postPoint->GetStepStatus() == fAtRestDoItProc ||
405  aTrack->GetTrackStatus() == fStopAndKill)
406  {
407  // save only hits with energy deposit (or -1 for geantino)
408  if (m_Hit->get_edep())
409  {
410  m_HitContainer->AddHit(layer_id, m_Hit);
411  if (m_SaveShower)
412  {
414  }
415  // ownership has been transferred to container, set to null
416  // so we will create a new hit for the next track
417  m_Hit = nullptr;
418  }
419  else
420  {
421  // if this hit has no energy deposit, just reset it for reuse
422  // this means we have to delete it in the dtor. If this was
423  // the last hit we processed the memory is still allocated
424  m_Hit->Reset();
425  }
426  }
427 
428  // return true to indicate the hit was used
429  return true;
430  }
431  else
432  {
433  return false;
434  }
435 
436  return false;
437 }
438 
439 //____________________________________________________________________________..
441 {
442  string hitnodename;
443  string absorbernodename;
444  if (m_Detector->SuperDetector() != "NONE")
445  {
446  hitnodename = "G4HIT_" + m_Detector->SuperDetector();
447  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->SuperDetector();
448  }
449  else
450  {
451  hitnodename = "G4HIT_" + m_Detector->GetName();
452  absorbernodename = "G4HIT_ABSORBER_" + m_Detector->GetName();
453  }
454 
455  //now look for the map and grab a pointer to it.
456  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
457  m_AbsorberhitContainer = findNode::getClass<PHG4HitContainer>(topNode, absorbernodename.c_str());
458 
459  // if we do not find the node it's messed up.
460  if (!m_HitContainer)
461  {
462  std::cout << "PHG4MvtxSteppingAction::SetTopNode - unable to find " << hitnodename << std::endl;
463  }
464  if (!m_AbsorberhitContainer)
465  {
466  if (Verbosity() > 0)
467  {
468  cout << "PHG4MvtxSteppingAction::SetTopNode - unable to find " << absorbernodename << endl;
469  }
470  }
471 }