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"))
76 G4TouchableHandle touch = aStep->GetPreStepPoint()->GetTouchableHandle();
100 if (whichactive == 2)
FindIndexZDC(touch, idx_j, idx_k);
101 if (whichactive == 1)
FindIndexSMD(touch, idx_j, idx_k);
104 G4double edep = aStep->GetTotalEnergyDeposit() /
GeV;
105 G4double eion = (aStep->GetTotalEnergyDeposit() - aStep->GetNonIonizingEnergyDeposit()) /
GeV;
106 G4double light_yield = 0;
109 const G4Track* aTrack = aStep->GetTrack();
114 edep = aTrack->GetKineticEnergy() /
GeV;
115 G4Track* killtrack =
const_cast<G4Track*
>(aTrack);
116 killtrack->SetTrackStatus(fStopAndKill);
123 bool geantino =
false;
124 if (aTrack->GetParticleDefinition()->GetPDGEncoding() == 0 &&
125 aTrack->GetParticleDefinition()->GetParticleName().find(
"geantino") != std::string::npos)
131 G4StepPoint* prePoint = aStep->GetPreStepPoint();
132 G4StepPoint* postPoint = aStep->GetPostStepPoint();
137 double charge = aTrack->GetParticleDefinition()->GetPDGCharge();
148 int pid = aTrack->GetParticleDefinition()->GetPDGEncoding();
150 const G4DynamicParticle* dypar = aTrack->GetDynamicParticle();
151 const G4ThreeVector& pdirect = dypar->GetMomentumDirection();
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)
162 G4double
E = dypar->GetTotalEnergy();
172 G4double
E = dypar->GetTotalEnergy();
173 G4double
P = dypar->GetTotalMomentum();
183 switch (prePoint->GetStepStatus())
198 m_Hit->
set_t(0, prePoint->GetGlobalTime() / nanosecond);
218 if (whichactive == -1)
228 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
242 if (whichactive == 1)
246 static bool once =
true;
248 if (once && edep > 0)
254 std::cout <<
"PHG4ZDCSteppingAction::UserSteppingAction::"
257 <<
" use scintillating light model at each Geant4 steps. "
260 << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetName()
262 <<
"Birk Constant = "
263 << aTrack->GetMaterialCutsCouple()->GetMaterial()->GetIonisation()->GetBirksConstant()
265 <<
"edep = " << edep <<
", "
268 <<
"light_yield = " << light_yield << std::endl;
288 if (postPoint->GetStepStatus() == fGeomBoundary ||
289 postPoint->GetStepStatus() == fWorldBoundary ||
290 postPoint->GetStepStatus() == fAtRestDoItProc ||
291 aTrack->GetTrackStatus() == fStopAndKill)
298 m_Hit->
set_t(1, postPoint->GetGlobalTime() / nanosecond);
321 if (G4VUserTrackInformation*
p = aTrack->GetUserInformation())
358 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_HitNodeName << std::endl;
362 if (!m_AbsorberHitContainer)
366 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_AbsorberNodeName << std::endl;
369 if (!m_SupportHitContainer)
373 std::cout <<
"PHG4ZDCSteppingAction::SetTopNode - unable to find " <<
m_SupportNodeName << std::endl;
385 else if (type ==
"G4HIT_ABSORBER")
390 else if (type ==
"G4HIT_SUPPORT")
395 std::cout <<
"Invalid output hit node type " << type << std::endl;
407 j = envelope->GetCopyNo();
408 k = (plate->GetCopyNo()) / 27;
420 int whichzdc = envelope->GetCopyNo();
422 j = scint->GetCopyNo() % 7;
423 k = scint->GetCopyNo() / 7;
425 if (whichzdc == 1) j += 7;
436 if (angle >= 90)
return 0;
437 for (
int i = 1;
i < 9;
i++)
441 std::array<double, 18> PMMAsub0 =
m_PMMA05[i - 1];
442 std::array<double, 18> PMMAsub1 =
m_PMMA05[
i];
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;
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;
463 if (E <
m_E[0])
return 0;
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);
475 for (
int i = 1;
i < 11;
i++)
479 std::array<double, 36> PMMAsub0 =
m_PMMA05E[i - 1];
480 std::array<double, 36> PMMAsub1 =
m_PMMA05E[
i];
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;
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;
496 std::cout <<
"out of range" << std::endl;