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