5 #include <phparameter/PHParameters.h>
17 #include <Geant4/G4DynamicParticle.hh>
18 #include <Geant4/G4IonisParamMat.hh>
19 #include <Geant4/G4Material.hh>
20 #include <Geant4/G4MaterialCutsCouple.hh>
21 #include <Geant4/G4ParticleDefinition.hh>
22 #include <Geant4/G4ReferenceCountedHandle.hh>
23 #include <Geant4/G4Step.hh>
24 #include <Geant4/G4StepPoint.hh>
25 #include <Geant4/G4StepStatus.hh>
26 #include <Geant4/G4String.hh>
27 #include <Geant4/G4SystemOfUnits.hh>
28 #include <Geant4/G4ThreeVector.hh>
29 #include <Geant4/G4TouchableHandle.hh>
30 #include <Geant4/G4Track.hh>
31 #include <Geant4/G4TrackStatus.hh>
32 #include <Geant4/G4Types.hh>
33 #include <Geant4/G4VPhysicalVolume.hh>
34 #include <Geant4/G4VTouchable.hh>
35 #include <Geant4/G4VUserTrackInformation.hh>
39 #include <gsl/gsl_randist.h>
40 #include <gsl/gsl_rng.h>
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"))
100 if (whichactive == 2)
FindIndexZDC(touch, idx_j, idx_k);
101 if (whichactive == 1)
FindIndexSMD(touch, idx_j, idx_k);
123 bool geantino =
false;
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)
216 if (whichactive == -1)
240 if (whichactive == 1)
244 static bool once =
true;
246 if (once && edep > 0)
252 std::cout <<
"PHG4ZDCSteppingAction::UserSteppingAction::"
255 <<
" use scintillating light model at each Geant4 steps. "
260 <<
"Birk Constant = "
263 <<
"edep = " << edep <<
", "
266 <<
"light_yield = " << light_yield << std::endl;
356 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_HitNodeName << std::endl;
360 if (!m_AbsorberHitContainer)
364 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_AbsorberNodeName << std::endl;
367 if (!m_SupportHitContainer)
371 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_SupportNodeName << std::endl;
383 else if (type ==
"G4HIT_ABSORBER")
388 else if (type ==
"G4HIT_SUPPORT")
393 std::cout <<
"Invalid output hit node type " << type << std::endl;
423 if (whichzdc == 1) j += 7;
434 if (angle >= 90)
return 0;
435 for (
int i = 1; i < 9; i++)
439 std::array<double, 18> PMMAsub0 =
m_PMMA05[i - 1];
440 std::array<double, 18> PMMAsub1 =
m_PMMA05[i];
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;
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;
461 if (E <
m_E[0])
return 0;
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);
473 for (
int i = 1; i < 11; i++)
477 std::array<double, 36> PMMAsub0 =
m_PMMA05E[i - 1];
478 std::array<double, 36> PMMAsub1 =
m_PMMA05E[i];
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;
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;
494 std::cout <<
"out of range" << std::endl;