ECCE @ EIC 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();
60  gsl_rng_set(RandomGenerator, seed);
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 //____________________________________________________________________________..
75 {
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  G4ThreeVector pdirect = dypar->GetMomentumDirection();
152  double dy = sqrt(2) / 2.;
153  double dz = sqrt(2) / 2.;
154  if (idx_j == 1) dz = -dz;
155  double CosTheta = pdirect.y() * dy + pdirect.z() * dz;
156  double angle = acos(CosTheta) * 180.0 / M_PI;
157  if (pid == 11 || pid == -11)
158  {
159  //find energy
160  G4double E = dypar->GetTotalEnergy();
161  //electron response here
162  double avg_ph = ZDCEResponse(E, angle);
163  avg_ph *= 0.16848;
164  //use Poisson Distribution here
165  int n_ph = gsl_ran_poisson(RandomGenerator, avg_ph);
166  light_yield += n_ph;
167  }
168  else
169  {
170  G4double E = dypar->GetTotalEnergy();
171  G4double P = dypar->GetTotalMomentum();
172  double beta = P / E;
173  double avg_ph = ZDCResponse(beta, angle);
174  avg_ph *= 0.16848;
175  int n_ph = gsl_ran_poisson(RandomGenerator, avg_ph);
176  light_yield += n_ph;
177  }
178  }
179  }
180  }
181  switch (prePoint->GetStepStatus())
182  {
183  case fGeomBoundary:
184  case fUndefined:
185  if (!m_Hit)
186  {
187  m_Hit = new PHG4Hitv1();
188  }
189 
190  /* Set hit location (space point) */
191  m_Hit->set_x(0, prePoint->GetPosition().x() / cm);
192  m_Hit->set_y(0, prePoint->GetPosition().y() / cm);
193  m_Hit->set_z(0, prePoint->GetPosition().z() / cm);
194 
195  /* Set hit time */
196  m_Hit->set_t(0, prePoint->GetGlobalTime() / nanosecond);
197 
198  //set the track ID
199  m_Hit->set_trkid(aTrack->GetTrackID());
200  /* set intial energy deposit */
201  m_Hit->set_edep(0);
202 
203  /* Now add the hit to the hit collection */
204  // here we do things which are different between scintillator and absorber hits
205  if (whichactive > 0)
206  {
208  m_Hit->set_eion(0);
209  m_Hit->set_light_yield(0); // for scintillator only, initialize light yields
210  /* Set hit location (tower index) */
211  m_Hit->set_index_k(idx_k);
212  m_Hit->set_index_j(idx_j);
213  }
214  else
215  {
216  if (whichactive == -1)
217  {
219  }
220  else
221  {
223  }
224  }
225  // here we set what is common for scintillator and absorber hits
227  {
228  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
229  {
230  m_Hit->set_trkid(pp->GetUserTrackId());
231  m_Hit->set_shower_id(pp->GetShower()->get_id());
232  m_CurrentShower = pp->GetShower();
233  }
234  }
235  break;
236  default:
237  break;
238  }
239 
240  if (whichactive == 1)
241  {
242  //light_yield = eion;
243  light_yield = GetVisibleEnergyDeposition(aStep); // for scintillator only, calculate light yields
244  static bool once = true;
245 
246  if (once && edep > 0)
247  {
248  once = false;
249 
250  if (Verbosity() > 0)
251  {
252  std::cout << "PHG4ZDCSteppingAction::UserSteppingAction::"
253  //
254  << m_Detector->GetName() << " - "
255  << " use scintillating light model at each Geant4 steps. "
256  << "First step: "
257  << "Material = "
258  << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
259  << ", "
260  << "Birk Constant = "
262  << ","
263  << "edep = " << edep << ", "
264  << "eion = " << eion
265  << ", "
266  << "light_yield = " << light_yield << std::endl;
267  }
268  }
269  }
270 
271  /* sum up the energy to get total deposited */
272 
274  if (whichactive > 0)
275  {
276  m_Hit->set_eion(m_Hit->get_eion() + eion);
277  m_Hit->set_light_yield(m_Hit->get_light_yield() + light_yield);
278  }
279 
280  // if any of these conditions is true this is the last step in
281  // this volume and we need to save the hit
282  // postPoint->GetStepStatus() == fGeomBoundary: track leaves this volume
283  // postPoint->GetStepStatus() == fWorldBoundary: track leaves this world
284  // (not sure if this will ever be the case)
285  // aTrack->GetTrackStatus() == fStopAndKill: track ends
286  if (postPoint->GetStepStatus() == fGeomBoundary ||
287  postPoint->GetStepStatus() == fWorldBoundary ||
288  postPoint->GetStepStatus() == fAtRestDoItProc ||
289  aTrack->GetTrackStatus() == fStopAndKill)
290  {
291  // Update exit values
292  m_Hit->set_x(1, postPoint->GetPosition().x() / cm);
293  m_Hit->set_y(1, postPoint->GetPosition().y() / cm);
294  m_Hit->set_z(1, postPoint->GetPosition().z() / cm);
295 
296  m_Hit->set_t(1, postPoint->GetGlobalTime() / nanosecond);
297 
298  // special case for geantinos
299  if (geantino)
300  {
301  m_Hit->set_edep(-1); // only energy=0 g4hits get dropped, this way geantinos survive the g4hit compression
302  if (whichactive > 0)
303  {
304  m_Hit->set_eion(-1);
305  m_Hit->set_light_yield(-1);
306  }
307  }
308  // save only hits with energy deposit (or -1 for geantino)
309  if (m_Hit->get_edep())
310  {
311  m_CurrentHitContainer->AddHit(layer_id, m_Hit);
312  if (m_CurrentShower)
313  {
315  }
316  // ownership has been transferred to container, set to null
317  if (m_Hit->get_edep() > 0 && (whichactive > 0 || absorbertruth > 0))
318  {
320  {
321  if (PHG4TrackUserInfoV1* pp = dynamic_cast<PHG4TrackUserInfoV1*>(p))
322  {
323  pp->SetKeep(1); // we want to keep the track
324  }
325  }
326  }
327  // so we will create a new hit for the next track
328  m_Hit = nullptr;
329  }
330  else
331  {
332  // if this hit has no energy deposit, just reset it for reuse
333  // this means we have to delete it in the dtor. If this was
334  // the last hit we processed the memory is still allocated
335  m_Hit->Reset();
336  }
337  }
338  return true;
339  }
340  else
341  {
342  return false;
343  }
344 }
345 
346 //____________________________________________________________________________..
348 {
349  //now look for the map and grab a pointer to it.
350  m_HitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_HitNodeName);
351  m_AbsorberHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_AbsorberNodeName);
352  m_SupportHitContainer = findNode::getClass<PHG4HitContainer>(topNode, m_SupportNodeName);
353  // if we do not find the node it's messed up.
354  if (!m_HitContainer)
355  {
356  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_HitNodeName << std::endl;
357  gSystem->Exit(1);
358  }
359  // this is perfectly fine if absorber hits are disabled
360  if (!m_AbsorberHitContainer)
361  {
362  if (Verbosity() > 0)
363  {
364  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_AbsorberNodeName << std::endl;
365  }
366  }
367  if (!m_SupportHitContainer)
368  {
369  if (Verbosity() > 0)
370  {
371  std::cout << "PHG4ZDCSteppingAction::SetTopNode - unable to find " << m_SupportNodeName << std::endl;
372  }
373  }
374 }
375 
376 void PHG4ZDCSteppingAction::SetHitNodeName(const std::string& type, const std::string& name)
377 {
378  if (type == "G4HIT")
379  {
381  return;
382  }
383  else if (type == "G4HIT_ABSORBER")
384  {
386  return;
387  }
388  else if (type == "G4HIT_SUPPORT")
389  {
391  return;
392  }
393  std::cout << "Invalid output hit node type " << type << std::endl;
394  gSystem->Exit(1);
395  return;
396 }
397 
398 //getting index using copyno
399 //if in ZDC PMMA fiber
401 {
402  G4VPhysicalVolume* envelope = touch->GetVolume(2); //Get the world
403  G4VPhysicalVolume* plate = touch->GetVolume(1); //Get the fiber plate
404 
405  j = envelope->GetCopyNo();
406  k = (plate->GetCopyNo()) / 27;
407 
408  return 0;
409 }
410 
412 {
413  int jshift = 10;
414  int kshift = 10;
415  G4VPhysicalVolume* envelope = touch->GetVolume(2); //Get the envelope
416  G4VPhysicalVolume* scint = touch->GetVolume(0); //Get the fiber plate
417 
418  int whichzdc = envelope->GetCopyNo();
419 
420  j = scint->GetCopyNo() % 7;
421  k = scint->GetCopyNo() / 7;
422 
423  if (whichzdc == 1) j += 7;
424  // shift the index to avoid conflict with the ZDC index
425  j += jshift;
426  k += kshift;
427 
428  return 0;
429 }
430 
431 double PHG4ZDCSteppingAction::ZDCResponse(double beta, double angle)
432 {
433  if (beta < m_BetaThersh) return 0;
434  if (angle >= 90) return 0;
435  for (int i = 1; i < 9; i++)
436  {
437  if (beta <= m_Beta[i])
438  {
439  std::array<double, 18> PMMAsub0 = m_PMMA05[i - 1];
440  std::array<double, 18> PMMAsub1 = m_PMMA05[i];
441  //find angle bin here and do 1D linear interpolation of angle
442  int Abin = (int) angle / 5;
443  if (Abin == 0) Abin = 1;
444  double avg_ph0 = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
445  double avg_ph1 = PMMAsub1[Abin - 1] + (PMMAsub1[Abin] - PMMAsub1[Abin - 1]) * (angle / 5 - Abin + 1);
446  if (avg_ph0 < 0) avg_ph0 = 0;
447  if (avg_ph1 < 0) avg_ph1 = 0;
448  //linear linear interpolation with beta
449  double avg_ph = avg_ph0 + (avg_ph1 - avg_ph0) * (beta - m_Beta[i - 1]) / (m_Beta[i] - m_Beta[i - 1]);
450  if (avg_ph < 0) avg_ph = 0;
451  //use poisson?
452  return avg_ph;
453  }
454  }
455 
456  return 0;
457 }
458 
460 {
461  if (E < m_E[0]) return 0;
462 
463  if (E >= 0.05)
464  {
465  std::array<double, 36> PMMAsub0 = m_PMMA05E[10];
466  int Abin = (int) angle / 5;
467  if (Abin == 0) Abin = 1;
468  double avg_ph = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
469  return avg_ph;
470  }
471  else
472  {
473  for (int i = 1; i < 11; i++)
474  {
475  if (E <= m_E[i])
476  {
477  std::array<double, 36> PMMAsub0 = m_PMMA05E[i - 1];
478  std::array<double, 36> PMMAsub1 = m_PMMA05E[i];
479 
480  int Abin = (int) angle / 5;
481  if (Abin == 0) Abin = 1;
482  double avg_ph0 = PMMAsub0[Abin - 1] + (PMMAsub0[Abin] - PMMAsub0[Abin - 1]) * (angle / 5 - Abin + 1);
483  double avg_ph1 = PMMAsub1[Abin - 1] + (PMMAsub1[Abin] - PMMAsub1[Abin - 1]) * (angle / 5 - Abin + 1);
484  if (avg_ph0 < 0) avg_ph0 = 0;
485  if (avg_ph1 < 0) avg_ph1 = 0;
486  //linear linear interpolation with E
487  double avg_ph = avg_ph0 + (avg_ph1 - avg_ph0) * (E - m_E[i - 1]) / (m_E[i] - m_E[i - 1]);
488  if (avg_ph < 0) avg_ph = 0;
489  //use poisson?
490  return avg_ph;
491  }
492  }
493  }
494  std::cout << "out of range" << std::endl;
495  return 0;
496 }