ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
EICG4LumiSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file EICG4LumiSteppingAction.cc
1 //____________________________________________________________________________..
2 //
3 // This is a working template for the Stepping Action which needs to be implemented
4 // for active detectors. Most of the code is error handling and access to the G4 objects
5 // and our data structures. It does not need any adjustment. The only thing you need to
6 // do is to add the properties of the G4Hits you want to save for later analysis
7 // This needs to be done in 2 places, G4Hits are generated when a G4 track enters a new
8 // volume (or is created). Here you give it an initial value. When the G4 track leaves
9 // the volume the final value needs to be set.
10 // The places to do this is marked by //implement your own here//
11 //
12 // As guidance you can look at the total (integrated over all steps in a volume) energy
13 // deposit which should always be saved.
14 // Additionally the total ionization energy is saved - this can be removed if you are not
15 // interested in this. Naturally you may want remove these comments in your version
16 //
17 //____________________________________________________________________________..
18 
20 
21 #include "EICG4LumiDetector.h"
22 #include "EICG4LumiSubsystem.h"
23 
24 #include <phparameter/PHParameters.h>
25 
27 
28 #include <g4main/PHG4Hit.h>
30 #include <g4main/PHG4Hitv1.h>
31 #include <g4main/PHG4Shower.h>
34 
35 #include <phool/getClass.h>
36 #include <phool/PHCompositeNode.h>
37 
38 #include <TSystem.h>
39 
40 #include <Geant4/G4NavigationHistory.hh>
41 #include <Geant4/G4ParticleDefinition.hh>
42 #include <Geant4/G4ReferenceCountedHandle.hh>
43 #include <Geant4/G4Step.hh>
44 #include <Geant4/G4StepPoint.hh>
45 #include <Geant4/G4StepStatus.hh>
46 #include <Geant4/G4String.hh>
47 #include <Geant4/G4SystemOfUnits.hh>
48 #include <Geant4/G4ThreeVector.hh>
49 #include <Geant4/G4TouchableHandle.hh>
50 #include <Geant4/G4Track.hh>
51 #include <Geant4/G4TrackStatus.hh>
52 #include <Geant4/G4Types.hh>
53 #include <Geant4/G4VPhysicalVolume.hh>
54 #include <Geant4/G4VTouchable.hh>
55 #include <Geant4/G4VUserTrackInformation.hh>
56 
57 #include <cmath>
58 #include <iostream>
59 #include <string>
60 #include <cstdlib>
61 #include <iomanip>
62 
63 class PHCompositeNode;
64 
65 //____________________________________________________________________________..
67  : PHG4SteppingAction(detector->GetName())
68  , m_Subsystem(subsys)
69  , m_Detector(detector)
70  , m_Params(parameters)
71  , m_HitContainerCAL(nullptr)
72  , m_HitContainerTracking(nullptr)
73  , m_HitContainerVirt(nullptr)
74  , m_Hit(nullptr)
75  , m_SaveShower(nullptr)
76  , m_SaveVolPre(nullptr)
77  , m_SaveVolPost(nullptr)
78  , m_SaveLightYieldFlag(m_Params->get_int_param("lightyield"))
79  , m_SaveTrackId(-1)
80  , m_SavePreStepStatus(-1)
81  , m_SavePostStepStatus(-1)
82  , m_ActiveFlag(m_Params->get_int_param("active"))
83  , m_BlackHoleFlag(m_Params->get_int_param("blackhole"))
84  , m_UseG4StepsFlag(m_Params->get_int_param("use_g4steps"))
85  , m_Tmin(m_Params->get_double_param("tmin") * ns)
86  , m_Tmax(m_Params->get_double_param("tmax") * ns)
87  , m_EdepSum(0)
88  , m_EabsSum(0)
89  , m_EionSum(0)
90 {
91 // G4 seems to have issues in the um range
92  m_Zmin -= copysign(m_Zmin, 1. / 1e6 * cm);
93  m_Zmax += copysign(m_Zmax, 1. / 1e6 * cm);
94 }
95 
96 //____________________________________________________________________________..
98 {
99  delete m_Hit;
100 }
101 
102 //____________________________________________________________________________..
103 // This is the implementation of the G4 UserSteppingAction
104 bool EICG4LumiSteppingAction::UserSteppingAction(const G4Step *aStep,bool was_used)
105 {
107  G4TouchableHandle touchpost = aStep->GetPostStepPoint()->GetTouchableHandle();
108  G4double light_yield = GetVisibleEnergyDeposition(aStep);
109 
110  // get volume of the current step
111  G4VPhysicalVolume *volume = touch->GetVolume();
112  G4LogicalVolume *volumeLogical = volume->GetLogicalVolume();
113  G4Material *mat = volumeLogical->GetMaterial();
114 
115  int activeMaterial = m_Detector->IsInDetector(volume);
116  int virtualMaterial = m_Detector->IsInVirtualDetector(volume);
117  int trackingMaterial = (mat->GetName().compareTo("G4_Si")==0) ? 1 : 0;
118 
119  if ( !activeMaterial && !virtualMaterial )
120  {
121  return false;
122  }
123  // collect energy and track length step by step
124  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
125  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
126  const G4Track *aTrack = aStep->GetTrack();
127  // if this detector stops everything, just put all kinetic energy into edep
128 
129  //G4ParticleDefinition *Pdef = aTrack->GetDefinition();
130  //double charge = Pdef->GetPDGCharge();
131  //if( abs(charge) != 1 ) { return false; }
132 
133  //G4cout<<"Name: "<<Pdef->GetParticleName()<<" Charge: "<<charge<<" P = "<<aTrack->GetMomentum()<<G4endl;
134 
135  if (m_BlackHoleFlag)
136  {
137  if ((!std::isfinite(m_Tmin) && !std::isfinite(m_Tmax)) ||
138  aTrack->GetGlobalTime() < m_Tmin ||
139  aTrack->GetGlobalTime() > m_Tmax)
140  {
141  edep = aTrack->GetKineticEnergy() / GeV;
142  G4Track *killtrack = const_cast<G4Track *>(aTrack);
143  killtrack->SetTrackStatus(fStopAndKill);
144  }
145 
146  }
147 
148  // int layer_id = m_Detector->get_Layer();
149  int layer_id = volume->GetCopyNo();
150  // To decompose into tower cells
151  //int idx_j = icopy >> 16;
152  //int idx_k = icopy & 0xFFFF;
153 
154  int layer_type = 1;
155  if (!m_ActiveFlag)
156  {
157  return false;
158  }
159  bool geantino = false;
160 
161  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
162  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") !=
163  std::string::npos) // this also accounts for "chargedgeantino"
164  {
165  geantino = true;
166  }
167  G4StepPoint *prePoint = aStep->GetPreStepPoint();
168  G4StepPoint *postPoint = aStep->GetPostStepPoint();
169 
170  switch (prePoint->GetStepStatus())
171  {
172  case fPostStepDoItProc:
174  {
175  // this is the okay case, fPostStepDoItProc called in a volume, not first thing inside
176  // a new volume, just proceed here
177  break;
178  }
179  else
180  {
181  // this is an impossible G4 Step print out diagnostic to help debug, not sure if
182  // this is still with us
183  std::cout << GetName() << ": New Hit for " << std::endl;
184  std::cout << "prestep status: "
186  << ", poststep status: "
188  << ", last pre step status: "
190  << ", last post step status: "
192  std::cout << "last track: " << m_SaveTrackId
193  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
194  std::cout << "phys pre vol: " << volume->GetName()
195  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
196  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
197  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
198  }
199  // These are the normal cases
200 
201  case fGeomBoundary:
202  case fUndefined:
203  if (!m_Hit)
204  {
205  m_Hit = new PHG4Hitv1();
206  }
207  m_Hit->set_layer((unsigned int) layer_id);
208  // here we set the entrance values in cm
209  //FindTowerIndexFromPosition(prePoint, idx_j, idx_k);
210  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
211  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
212  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
213 
214  m_Hit->set_px(0, prePoint->GetMomentum().x() / GeV);
215  m_Hit->set_py(0, prePoint->GetMomentum().y() / GeV);
216  m_Hit->set_pz(0, prePoint->GetMomentum().z() / GeV);
217 
218  // time in ns
219  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
220 
221  // set the tower index
222  //m_Hit->set_index_j(idx_j);
223  //m_Hit->set_index_k(idx_k);
224 
225  // set the track ID
226  m_Hit->set_trkid(aTrack->GetTrackID());
227  m_SaveTrackId = aTrack->GetTrackID();
228  // set the initial energy deposit
229  m_EdepSum = 0;
230  m_EabsSum = 0;
231 
232  m_Hit->set_edep(0);
233  if (!geantino && !m_BlackHoleFlag)
234  {
235  m_Hit->set_eion(0);
236  }
237 
239  {
241  }
242 
244  {
245  if (PHG4TrackUserInfoV1 *pp = dynamic_cast<PHG4TrackUserInfoV1 *>(p))
246  {
247  m_Hit->set_trkid(pp->GetUserTrackId());
248  m_Hit->set_shower_id(pp->GetShower()->get_id());
249  m_SaveShower = pp->GetShower();
250  }
251  }
252  if (!hasMotherSubsystem() && (m_Hit->get_z(0) * cm > m_Zmax || m_Hit->get_z(0) * cm < m_Zmin))
253  {
254  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
255  << " PHG4CylinderSteppingAction: Entry hit z " << m_Hit->get_z(0) * cm
256  << " outside acceptance, zmin " << m_Zmin
257  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
258  }
259  break;
260  default:
261  break;
262  }
263 
264  if (!m_Hit || !std::isfinite(m_Hit->get_x(0)))
265  {
266  std::cout << GetName() << ": hit was not created" << std::endl;
267  std::cout << "prestep status: "
269  << ", poststep status: "
271  << ", last pre step status: "
273  << ", last post step status: "
275  std::cout << "last track: " << m_SaveTrackId
276  << ", current trackid: " << aTrack->GetTrackID() << std::endl;
277  std::cout << "phys pre vol: " << volume->GetName()
278  << " post vol : " << touchpost->GetVolume()->GetName() << std::endl;
279  std::cout << " previous phys pre vol: " << m_SaveVolPre->GetName()
280  << " previous phys post vol: " << m_SaveVolPost->GetName() << std::endl;
281  // This is fatal - a hit from nowhere. This needs to be looked at and fixed
282  gSystem->Exit(1);
283  }
284  m_SavePostStepStatus = postPoint->GetStepStatus();
285  // check if track id matches the initial one when the hit was created
286  if (aTrack->GetTrackID() != m_SaveTrackId)
287  {
288  std::cout << GetName() << ": hits do not belong to the same track" << std::endl;
289  std::cout << "saved track: " << m_SaveTrackId
290  << ", current trackid: " << aTrack->GetTrackID()
291  << ", prestep status: " << prePoint->GetStepStatus()
292  << ", previous post step status: " << m_SavePostStepStatus << std::endl;
293  gSystem->Exit(1);
294  }
295 
296  // We need to cache a few things from one step to the next
297  // to identify impossible hits and subsequent debugging printout
298  m_SavePreStepStatus = prePoint->GetStepStatus();
299  m_SavePostStepStatus = postPoint->GetStepStatus();
301  m_SaveVolPost = touchpost->GetVolume();
302 
303  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
304  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
305  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
306  m_Hit->set_px(1, postPoint->GetMomentum().x() / GeV);
307  m_Hit->set_py(1, postPoint->GetMomentum().y() / GeV);
308  m_Hit->set_pz(1, postPoint->GetMomentum().z() / GeV);
309 
310  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
311  //sum up the energy to get total deposited
313  m_EdepSum += edep;
314  if (!hasMotherSubsystem() && (m_Hit->get_z(1) * cm > m_Zmax || m_Hit->get_z(1) * cm < m_Zmin))
315  {
316  std::cout << m_Detector->SuperDetector() << std::setprecision(9)
317  << " PHG4CylinderSteppingAction: Exit hit z " << m_Hit->get_z(1) * cm
318  << " outside acceptance zmin " << m_Zmin
319  << ", zmax " << m_Zmax << ", layer: " << layer_id << std::endl;
320  }
321  if (geantino)
322  {
323  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
324  m_Hit->set_eion(-1);
325  m_Hit->set_light_yield(-1);
326 
327  }
328  else
329  {
330  if (!m_BlackHoleFlag)
331  {
332  m_Hit->set_eion(m_Hit->get_eion() + eion);
334  m_Hit->set_light_yield(m_Hit->get_light_yield()+light_yield);
335  }
336  }
337 
338  if (edep > 0 || m_SaveAllHitsFlag)
339  {
341  {
342  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
343  {
344  pp->SetKeep(1); // we want to keep the track
345  }
346  }
347  }
348 
349  if (postPoint->GetStepStatus() == fGeomBoundary ||
350  postPoint->GetStepStatus() == fWorldBoundary ||
351  postPoint->GetStepStatus() == fAtRestDoItProc ||
352  aTrack->GetTrackStatus() == fStopAndKill ||
353  m_UseG4StepsFlag > 0)
354  {
355  // save only hits with energy deposit (or geantino)
356  if (m_Hit->get_edep() || m_SaveAllHitsFlag)
357  {
358  m_Hit->set_layer(layer_id);
359  m_Hit->set_hit_type(layer_type);
360  if( activeMaterial )
361  {
362  if( trackingMaterial ) {
363  m_HitContainerTracking->AddHit(layer_id, m_Hit);
364  }
365  else {
366  m_HitContainerCAL->AddHit(layer_id, m_Hit);
367  }
368  }
369  if( virtualMaterial )
370  {
371  m_HitContainerVirt->AddHit(layer_id, m_Hit);
372  }
373 
374  if (m_SaveShower)
375  {
377  }
378  m_Hit = nullptr;
379  }
380  else
381  {
382  m_Hit->Reset();
383  }
384  }
385 
386  // return true to indicate the hit was used
387  return true;
388 }
389 
390 //____________________________________________________________________________..
392 {
393  //std::string hitnodename = "G4HIT_" + m_Detector->GetName();
394 
395  // now look for the map and grab a pointer to it.
396  m_HitContainerCAL = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeNameCAL);
397  m_HitContainerTracking = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeNameTracking);
398  m_HitContainerVirt = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeNameVirt);
399 
400  // if we do not find the node we need to make it.
401  if (!m_HitContainerCAL)
402  {
403  std::cout << PHWHERE << " EICG4LumiSteppingAction::SetTopNode - unable to find "
404  << m_HitNodeNameCAL << std::endl;
405 
406  gSystem->Exit(1);
407  }
408  if (!m_HitContainerTracking)
409  {
410  std::cout << PHWHERE << " EICG4LumiSteppingAction::SetTopNode - unable to find "
411  << m_HitNodeNameTracking << std::endl;
412 
413  gSystem->Exit(1);
414  }
415  if (!m_HitContainerVirt)
416  {
417  std::cout << PHWHERE << " EICG4LumiSteppingAction::SetTopNode - unable to find "
418  << m_HitNodeNameVirt << std::endl;
419 
420  gSystem->Exit(1);
421  }
422 }
423 
424 //__________________________________________________________________________..
426 {
428  {
429  return true;
430  }
431  return false;
432 }