3 #include <phparameter/PHParameterInterface.h>
29 #include <gsl/gsl_const_mksa.h>
53 static constexpr
double cm = 1.0;
57 static constexpr
double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1
e-7;
60 static constexpr
double maxHitLength = 1. *
cm;
63 static constexpr
double halflength_tpc = 105.5 *
cm;
66 static constexpr
double begin_CM = 20. *
cm;
67 static constexpr
double end_CM = 78. *
cm;
70 static constexpr
double halfwidth_CM = 0.5 *
cm;
73 std::optional<TVector3> central_membrane_intersection(TVector3
start, TVector3 direction)
75 const double end = start.z() > 0 ? halfwidth_CM : -halfwidth_CM;
76 const double dist = end - start.z();
79 if (direction.z() == 0)
return std::nullopt;
82 if (dist * direction.z() < 0)
return std::nullopt;
84 const double direction_scale = dist / direction.z();
85 return start + direction * direction_scale;
89 std::optional<TVector3> endcap_intersection(TVector3 start, TVector3 direction)
91 const double end = start.z() > 0 ? halflength_tpc : -halflength_tpc;
92 const double dist = end - start.z();
95 if (direction.z() == 0)
return std::nullopt;
98 if (dist * direction.z() < 0)
return std::nullopt;
100 const double direction_scale = dist / direction.z();
101 return start + direction * direction_scale;
105 std::optional<TVector3> cylinder_line_intersection(TVector3
s, TVector3
v,
double radius)
107 const double R2 =
square(radius);
112 const double b = 2 * (v.x() * s.x() + v.y() * s.y());
115 const double rootterm =
square(b) - 4 * a *
c;
122 if (rootterm < 0 || a == 0)
return std::nullopt;
125 const double sqrtterm = std::sqrt(rootterm);
126 const double t1 = (-b + sqrtterm) / (2 * a);
127 const double t2 = (-b - sqrtterm) / (2 * a);
133 const double& min_t = (t2 < t1 && t2 > 0) ? t2 : t1;
134 return s + v * min_t;
138 std::optional<TVector3> field_cage_intersection(TVector3 start, TVector3 direction)
140 const auto ofc_strike = cylinder_line_intersection(start, direction, end_CM);
141 const auto ifc_strike = cylinder_line_intersection(start, direction, begin_CM);
144 if (!ifc_strike)
return ofc_strike;
145 if (!ofc_strike)
return ifc_strike;
148 const auto ifc_dist = (ifc_strike->Z() - start.Z()) / direction.Z();
149 const auto ofc_dist = (ofc_strike->Z() - start.Z()) / direction.Z();
152 return (ofc_dist > 0) ? ofc_strike : std::nullopt;
153 else if (ofc_dist < 0)
156 return (ifc_dist < ofc_dist) ? ifc_strike : ofc_strike;
160 inline std::ostream&
operator<<(std::ostream& out,
const TVector3& vector)
162 out <<
"( " << vector.x() <<
", " << vector.y() <<
", " << vector.z() <<
")";
180 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
183 std::cout <<
"Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
189 std::cout <<
PHWHERE <<
"DST Node missing, aborting." << std::endl;
199 auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode,
hitnodename.c_str());
202 std::cout <<
Name() <<
" Could not locate G4HIT node " <<
hitnodename << std::endl;
216 std::cout <<
PHWHERE <<
"DST Node missing, aborting." << std::endl;
240 std::cout <<
"PHG4TpcDirectLaser::InitRun - phi steps: " <<
nPhiSteps <<
" min: " <<
minPhi <<
" max: " <<
maxPhi << std::endl;
241 std::cout <<
"PHG4TpcDirectLaser::InitRun - theta steps: " <<
nThetaSteps <<
" min: " <<
minTheta <<
" max: " <<
maxTheta << std::endl;
242 std::cout <<
"PHG4TpcDirectLaser::InitRun - nTotalSteps: " <<
nTotalSteps << std::endl;
244 std::cout <<
"PHG4TpcDirectLaser::InitRun - electrons_per_cm: " <<
electrons_per_cm << std::endl;
245 std::cout <<
"PHG4TpcDirectLaser::InitRun - electrons_per_gev " <<
electrons_per_gev << std::endl;
254 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
287 static constexpr
double Ne_dEdx = 1.56;
288 static constexpr
double CF4_dEdx = 7.00;
289 static constexpr
double Ne_NTotal = 43;
290 static constexpr
double CF4_NTotal = 100;
291 static constexpr
double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
292 static constexpr
double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
293 static constexpr
double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
305 if (n < 0 || max < min)
307 std::cout <<
PHWHERE <<
" - invalid" << std::endl;
319 if (n < 0 || max < min)
321 std::cout <<
PHWHERE <<
" - invalid" << std::endl;
339 const TVector3 position_base(60 * cm, 0., halflength_tpc);
342 for (
int i = 0; i < 8; ++i)
387 std::cout <<
"PHG4TpcDirectLaser::AimToThetaPhi - theta: " << theta <<
" phi: " << phi << std::endl;
405 std::cout <<
"PHG4TpcDirectLaser::AimToPatternStep - step: " << n <<
"/" <<
nTotalSteps << std::endl;
430 std::cout <<
PHWHERE <<
"invalid g4hit container. aborting" << std::endl;
439 TVector3
dir(0, 0, direction);
442 dir.RotateY(theta * direction);
446 dir.RotateZ(laser.
m_phi);
451 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack - position: " <<
pos <<
" direction: " << dir << std::endl;
455 static constexpr
double total_momentum = 1;
465 const auto vertex =
new PHG4VtxPoint_t(
pos.x(),
pos.y(),
pos.z(), 0, vtxid);
472 auto particle =
new PHG4Particle_t();
477 particle->set_px(total_momentum * dir.x());
478 particle->set_py(total_momentum * dir.y());
479 particle->set_pz(total_momentum * dir.z());
493 track.
set_px(total_momentum * dir.x());
494 track.
set_py(total_momentum * dir.y());
495 track.
set_pz(total_momentum * dir.z());
502 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack - position: " <<
pos <<
" direction: " << dir << std::endl;
512 const auto plane_strike = (
pos.z() * dir.z() > 0) ? endcap_intersection(
pos, dir) : central_membrane_intersection(
pos, dir);
515 const auto fc_strike = field_cage_intersection(
pos, dir);
518 if (!(plane_strike || fc_strike))
return;
522 const TVector3& strike = (fc_strike && (!plane_strike || fc_strike->z() / dir.z() < plane_strike->z() / dir.z())) ? *fc_strike : *plane_strike;
526 double fullLength = delta.Mag();
527 int nHitSteps = fullLength / maxHitLength + 1;
529 TVector3 start =
pos;
530 TVector3 end =
start;
531 TVector3
step = dir * (maxHitLength / (dir.Mag()));
532 double stepLength = 0;
536 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack -"
537 <<
" fullLength: " << fullLength
538 <<
" nHitSteps: " << nHitSteps
542 for (
int i = 0; i < nHitSteps; i++)
545 if (i + 1 == nHitSteps)
550 stepLength = delta.Mag();
556 stepLength = step.Mag();
560 auto hit =
new PHG4Hit_t;
561 hit->set_trkid(trackid);
565 hit->set_x(0, start.X() /
cm);
566 hit->set_y(0, start.Y() /
cm);
567 hit->set_z(0, start.Z() /
cm);
568 hit->set_t(0, (start -
pos).Mag() / speed_of_light);
570 hit->set_x(1, end.X() /
cm);
571 hit->set_y(1, end.Y() /
cm);
572 hit->set_z(1, end.Z() /
cm);
573 hit->set_t(1, (end -
pos).Mag() / speed_of_light);
576 hit->set_px(0, dir.X());
577 hit->set_py(0, dir.Y());
578 hit->set_pz(0, dir.Z());
580 hit->set_px(1, dir.X());
581 hit->set_py(1, dir.Y());
582 hit->set_pz(1, dir.Z());
586 hit->set_eion(totalE);
587 hit->set_edep(totalE);