Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ZDCSteppingAction.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ZDCSteppingAction.cc
2 
3 #include "PHG4ZDCDetector.h"
4 
5 #include <phparameter/PHParameters.h>
6 
7 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4Hitv1.h>
10 #include <g4main/PHG4Shower.h>
11 #include <g4main/PHG4SteppingAction.h> // for PHG4SteppingAction
13 
14 #include <phool/PHRandomSeed.h>
15 #include <phool/getClass.h>
16 
17 #include <Geant4/G4DynamicParticle.hh> // for G4DynamicParticle
18 #include <Geant4/G4IonisParamMat.hh> // for G4IonisParamMat
19 #include <Geant4/G4Material.hh> // for G4Material
20 #include <Geant4/G4MaterialCutsCouple.hh>
21 #include <Geant4/G4ParticleDefinition.hh> // for G4ParticleDefinition
22 #include <Geant4/G4ReferenceCountedHandle.hh> // for G4ReferenceCountedHandle
23 #include <Geant4/G4Step.hh>
24 #include <Geant4/G4StepPoint.hh> // for G4StepPoint
25 #include <Geant4/G4StepStatus.hh> // for fGeomBoundary, fAtRest...
26 #include <Geant4/G4String.hh> // for G4String
27 #include <Geant4/G4SystemOfUnits.hh>
28 #include <Geant4/G4ThreeVector.hh> // for G4ThreeVector
29 #include <Geant4/G4TouchableHandle.hh>
30 #include <Geant4/G4Track.hh> // for G4Track
31 #include <Geant4/G4TrackStatus.hh> // for fStopAndKill
32 #include <Geant4/G4Types.hh> // for G4double
33 #include <Geant4/G4VPhysicalVolume.hh> // for G4VPhysicalVolume
34 #include <Geant4/G4VTouchable.hh> // for G4VTouchable
35 #include <Geant4/G4VUserTrackInformation.hh> // for G4VUserTrackInformation
36 
37 #include <TSystem.h>
38 
39 #include <gsl/gsl_randist.h>
40 #include <gsl/gsl_rng.h>
41 
42 #include <array>
43 #include <cmath>
44 #include <iostream>
45 #include <string> // for basic_string, operator+
46 
47 class PHCompositeNode;
48 
49 //____________________________________________________________________________..
51  : PHG4SteppingAction(detector->GetName())
52  , m_Detector(detector)
53  , m_Params(parameters)
54  , m_IsActiveFlag(m_Params->get_int_param("active"))
55  , absorbertruth(m_Params->get_int_param("absorberactive"))
56  , m_IsBlackHole(m_Params->get_int_param("blackhole"))
57 {
58  RandomGenerator = gsl_rng_alloc(gsl_rng_mt19937);
59  unsigned int seed = PHRandomSeed();
61 }
62 
64 {
65  // if the last hit was a zero energie deposit hit, it is just reset
66  // and the memory is still allocated, so we need to delete it here
67  // if the last hit was saved, hit is a nullptr pointer which are
68  // legal to delete (it results in a no operation)
69  delete m_Hit;
70  gsl_rng_free(RandomGenerator);
71 }
72 
73 //____________________________________________________________________________..
74 bool PHG4ZDCSteppingAction::UserSteppingAction(const G4Step* aStep, bool)
75 {
76  G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
77  G4VPhysicalVolume* volume = touch->GetVolume();
78 
79  // m_Detector->IsInZDC(volume)
80  // returns
81  // 0 is outside of ZDC
82  // 1 is inside scintillator
83  // -1 is inside absorber or support structure (dead material)
84 
85  int whichactive = m_Detector->IsInZDC(volume);
86 
87  if (!whichactive)
88  {
89  return false;
90  }
91 
92  int layer_id = m_Detector->get_Layer();
93  int idx_k = -1;
94  int idx_j = -1;
95 
96  if (whichactive > 0) // in scintillator or fiber
97  {
98  /* Find indices of scintillator / tower containing this step */
99  //FindIndex(touch, idx_j, idx_k);
100  if (whichactive == 2) FindIndexZDC(touch, idx_j, idx_k);
101  if (whichactive == 1) FindIndexSMD(touch, idx_j, idx_k);
102  }
103  /* Get energy deposited by this step */
104  G4double edep = aStep->GetTotalEnergyDeposit() / GeV;
105  G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) / GeV;
106  G4double light_yield = 0;
107 
108  /* Get pointer to associated Geant4 track */
109  const G4Track* aTrack = aStep->GetTrack();
110 
111  // if this block stops everything, just put all kinetic energy into edep
112  if (m_IsBlackHole)
113  {
114  edep = aTrack->GetKineticEnergy() / GeV;
115  G4Track* killtrack = const_cast<G4Track*>(aTrack);
116  killtrack->SetTrackStatus(fStopAndKill);
117  }
118 
119  /* Make sure we are in a volume */
120  if (m_IsActiveFlag)
121  {
122  /* Check if particle is 'geantino' */
123  bool geantino = false;
124  if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
125  aTrack->GetParticleDefinition()->GetParticleName().find("geantino") != std::string::npos)
126  {
127  geantino = true;
128  }
129 
130  /* Get Geant4 pre- and post-step points */
131  G4StepPoint* prePoint = aStep->GetPreStepPoint();
132  G4StepPoint* postPoint = aStep->GetPostStepPoint();
133 
134  //if prepoint is in fiber
135  if (whichactive > 1)
136  {
137  double charge = aTrack->GetParticleDefinition()->GetPDGCharge();
138  //if charged particle
139  if (charge != 0)
140  {
141  //check if prepoint in active volume & postpoint out of active volume
142  G4VPhysicalVolume* postvolume = postPoint->GetTouchableHandle()->GetVolume();
143  int postactive = m_Detector->IsInZDC(postvolume);
144  //postpoint outside fiber
145  if (!postactive)
146  {
147  //get particle information here
148  int pid = aTrack->GetParticleDefinition()->GetPDGEncoding();
149  //calculate incidence angle
150  const G4DynamicParticle* dypar = aTrack->GetDynamicParticle();
151  const G4ThreeVector& pdirect = dypar->GetMomentumDirection();
152  // this triggers cppcheck, the code is good and the warning is suppressed
153  // cppcheck-suppress [duplicateAssignExpression, unmatchedSuppression]
154  double dy = sqrt(2) / 2.;
155  double dz = sqrt(2) / 2.;
156  if (idx_j == 1) dz = -dz;
157  double CosTheta = pdirect.y() * dy + pdirect.z() * dz;
158  double angle = acos(CosTheta) * 180.0 / M_PI;
159  if (pid == 11 || pid == -11)
160  {
161  //find energy
162  G4double E = dypar->GetTotalEnergy();
163  //electron response here
164  double avg_ph = ZDCEResponse(E, angle);
165  avg_ph *= 0.16848;
166  //use Poisson Distribution here
167  int n_ph = gsl_ran_poisson(RandomGenerator, avg_ph);
168  light_yield += n_ph;
169  }
170  else
171  {
172  G4double E = dypar->GetTotalEnergy();
173  G4double P = dypar->GetTotalMomentum();
174  double beta = P / E;
175  double avg_ph = ZDCResponse(beta, angle);
176  avg_ph *= 0.16848;
177  int n_ph = gsl_ran_poisson(RandomGenerator, avg_ph);
178  light_yield += n_ph;
179  }
180  }
181  }
182  }
183  switch (prePoint->GetStepStatus())
184  {
185  case fGeomBoundary:
186  case fUndefined:
187  if (!m_Hit)
188  {
189  m_Hit = new PHG4Hitv1();
190  }
191 
192  /* Set hit location (space point) */
193  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
194  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
195  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
196 
197  /* Set hit time */
198  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
199 
200  //set the track ID
201  m_Hit->set_trkid(aTrack->GetTrackID());
202  /* set intial energy deposit */
203  m_Hit->set_edep(0);
204 
205  /* Now add the hit to the hit collection */
206  // here we do things which are different between scintillator and absorber hits
207  if (whichactive > 0)
208  {
210  m_Hit->set_eion(0);
211  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
212  /* Set hit location (tower index) */
213  m_Hit->set_index_k(idx_k);
214  m_Hit->set_index_j(idx_j);
215  }
216  else
217  {
218  if (whichactive == -1)
219  {
221  }
222  else
223  {
225  }
226  }
227  // here we set what is common for scintillator and absorber hits
228  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
229  {
230  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
231  {
232  m_Hit->set_trkid(pp->GetUserTrackId());
233  m_Hit->set_shower_id(pp->GetShower()->get_id());
234  m_CurrentShower = pp->GetShower();
235  }
236  }
237  break;
238  default:
239  break;
240  }
241 
242  if (whichactive == 1)
243  {
244  //light_yield = eion;
245  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
246  static bool once = true;
247 
248  if (once && edep > 0)
249  {
250  once = false;
251 
252  if (Verbosity() > 0)
253  {
254  std::cout << "PHG4ZDCSteppingAction::UserSteppingAction::"
255  //
256  << m_Detector->GetName() << " - "
257  << " use scintillating light model at each Geant4 steps. "
258  << "First step: "
259  << "Material = "
260  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
261  << ", "
262  << "Birk Constant = "
263  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetIonisation()->GetBirksConstant()
264  << ","
265  << "edep = " << edep << ", "
266  << "eion = " << eion
267  << ", "
268  << "light_yield = " << light_yield << std::endl;
269  }
270  }
271  }
272 
273  /* sum up the energy to get total deposited */
274 
275  m_Hit->set_edep(m_Hit->get_edep() + edep);
276  if (whichactive > 0)
277  {
278  m_Hit->set_eion(m_Hit->get_eion() + eion);
279  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
280  }
281 
282  // if any of these conditions is true this is the last step in
283  // this volume and we need to save the hit
284  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
285  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
286  // (not sure if this will ever be the case)
287  // aTrack->GetTrackStatus() == fStopAndKill: track ends
288  if (postPoint->GetStepStatus() == fGeomBoundary ||
289  postPoint->GetStepStatus() == fWorldBoundary ||
290  postPoint->GetStepStatus() == fAtRestDoItProc ||
291  aTrack->GetTrackStatus() == fStopAndKill)
292  {
293  // Update exit values
294  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
295  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
296  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
297 
298  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
299 
300  // special case for geantinos
301  if (geantino)
302  {
303  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
304  if (whichactive > 0)
305  {
306  m_Hit->set_eion(-1);
307  m_Hit->set_light_yield(-1);
308  }
309  }
310  // save only hits with energy deposit (or -1 for geantino)
311  if (m_Hit->get_edep())
312  {
313  m_CurrentHitContainer->AddHit(layer_id, m_Hit);
314  if (m_CurrentShower)
315  {
317  }
318  // ownership has been transferred to container, set to null
319  if (m_Hit->get_edep() > 0 && (whichactive > 0 || absorbertruth > 0))
320  {
321  if (G4VUserTrackInformation* p = aTrack->GetUserInformation())
322  {
323  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
324  {
325  pp->SetKeep(1); // we want to keep the track
326  }
327  }
328  }
329  // so we will create a new hit for the next track
330  m_Hit = nullptr;
331  }
332  else
333  {
334  // if this hit has no energy deposit, just reset it for reuse
335  // this means we have to delete it in the dtor. If this was
336  // the last hit we processed the memory is still allocated
337  m_Hit->Reset();
338  }
339  }
340  return true;
341  }
342  else
343  {
344  return false;
345  }
346 }
347 
348 //____________________________________________________________________________..
350 {
351  //now look for the map and grab a pointer to it.
352  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
353  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
354  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
355  // if we do not find the node it's messed up.
356  if (!m_HitContainer)
357  {
358  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
359  gSystem->Exit(1);
360  }
361  // this is perfectly fine if absorber hits are disabled
362  if (!m_AbsorberHitContainer)
363  {
364  if (Verbosity() > 0)
365  {
366  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
367  }
368  }
369  if (!m_SupportHitContainer)
370  {
371  if (Verbosity() > 0)
372  {
373  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
374  }
375  }
376 }
377 
379 {
380  if (type == "G4HIT")
381  {
383  return;
384  }
385  else if (type == "G4HIT_ABSORBER")
386  {
388  return;
389  }
390  else if (type == "G4HIT_SUPPORT")
391  {
393  return;
394  }
395  std::cout << "Invalid output hit node type " << type << std::endl;
396  gSystem->Exit(1);
397  return;
398 }
399 
400 //getting index using copyno
401 //if in ZDC PMMA fiber
402 int PHG4ZDCSteppingAction::FindIndexZDC(G4TouchableHandle& touch, int& j, int& k)
403 {
404  G4VPhysicalVolume* envelope = touch->GetVolume(2); //Get the world
405  G4VPhysicalVolume* plate = touch->GetVolume(1); //Get the fiber plate
406 
407  j = envelope->GetCopyNo();
408  k = (plate->GetCopyNo()) / 27;
409 
410  return 0;
411 }
412 
413 int PHG4ZDCSteppingAction::FindIndexSMD(G4TouchableHandle& touch, int& j, int& k)
414 {
415  int jshift = 10;
416  int kshift = 10;
417  G4VPhysicalVolume* envelope = touch->GetVolume(2); //Get the envelope
418  G4VPhysicalVolume* scint = touch->GetVolume(0); //Get the fiber plate
419 
420  int whichzdc = envelope->GetCopyNo();
421 
422  j = scint->GetCopyNo() % 7;
423  k = scint->GetCopyNo() / 7;
424 
425  if (whichzdc == 1) j += 7;
426  // shift the index to avoid conflict with the ZDC index
427  j += jshift;
428  k += kshift;
429 
430  return 0;
431 }
432 
433 double PHG4ZDCSteppingAction::ZDCResponse(double beta, double angle)
434 {
435  if (beta < m_BetaThersh) return 0;
436  if (angle >= 90) return 0;
437  for (int i = 1; i < 9; i++)
438  {
439  if (beta <= m_Beta[i])
440  {
441  std::array<double, 18> PMMAsub0 = m_PMMA05[i - 1];
442  std::array<double, 18> PMMAsub1 = m_PMMA05[i];
443  //find angle bin here and do 1D linear interpolation of angle
444  int Abin = (int) angle / 5;
445  if (Abin == 0) Abin = 1;
446  double avg_ph0 = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
447  double avg_ph1 = PMMAsub1[Abin - 1] + (PMMAsub1[Abin] - PMMAsub1[Abin - 1]) * (angle / 5 - Abin + 1);
448  if (avg_ph0 < 0) avg_ph0 = 0;
449  if (avg_ph1 < 0) avg_ph1 = 0;
450  //linear linear interpolation with beta
451  double avg_ph = avg_ph0 + (avg_ph1 - avg_ph0) * (beta - m_Beta[i - 1]) / (m_Beta[i] - m_Beta[i - 1]);
452  if (avg_ph < 0) avg_ph = 0;
453  //use poisson?
454  return avg_ph;
455  }
456  }
457 
458  return 0;
459 }
460 
461 double PHG4ZDCSteppingAction::ZDCEResponse(double E, double angle)
462 {
463  if (E < m_E[0]) return 0;
464 
465  if (E >= 0.05)
466  {
467  std::array<double, 36> PMMAsub0 = m_PMMA05E[10];
468  int Abin = (int) angle / 5;
469  if (Abin == 0) Abin = 1;
470  double avg_ph = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
471  return avg_ph;
472  }
473  else
474  {
475  for (int i = 1; i < 11; i++)
476  {
477  if (E <= m_E[i])
478  {
479  std::array<double, 36> PMMAsub0 = m_PMMA05E[i - 1];
480  std::array<double, 36> PMMAsub1 = m_PMMA05E[i];
481 
482  int Abin = (int) angle / 5;
483  if (Abin == 0) Abin = 1;
484  double avg_ph0 = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
485  double avg_ph1 = PMMAsub1[Abin - 1] + (PMMAsub1[Abin] - PMMAsub1[Abin - 1]) * (angle / 5 - Abin + 1);
486  if (avg_ph0 < 0) avg_ph0 = 0;
487  if (avg_ph1 < 0) avg_ph1 = 0;
488  //linear linear interpolation with E
489  double avg_ph = avg_ph0 + (avg_ph1 - avg_ph0) * (E - m_E[i - 1]) / (m_E[i] - m_E[i - 1]);
490  if (avg_ph < 0) avg_ph = 0;
491  //use poisson?
492  return avg_ph;
493  }
494  }
495  }
496  std::cout << "out of range" << std::endl;
497  return 0;
498 }