3 #include <phparameter/PHParameterInterface.h>
31 #include <gsl/gsl_const_mksa.h>
55 static constexpr
double cm = 1.0;
59 static constexpr
double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1
e-7;
62 static constexpr
double maxHitLength = 1. *
cm;
65 static constexpr
double halflength_tpc = 105.5 *
cm;
68 static constexpr
double begin_CM = 20. *
cm;
69 static constexpr
double end_CM = 78. *
cm;
72 static constexpr
double halfwidth_CM = 0.5 *
cm;
75 std::optional<TVector3> central_membrane_intersection(TVector3
start, TVector3 direction)
77 const double end = start.z() > 0 ? halfwidth_CM : -halfwidth_CM;
78 const double dist = end - start.z();
81 if (direction.z() == 0)
return std::nullopt;
84 if (dist * direction.z() < 0)
return std::nullopt;
86 const double direction_scale = dist / direction.z();
87 return start + direction * direction_scale;
91 std::optional<TVector3> endcap_intersection(TVector3 start, TVector3 direction)
93 const double end = start.z() > 0 ? halflength_tpc : -halflength_tpc;
94 const double dist = end - start.z();
97 if (direction.z() == 0)
return std::nullopt;
100 if (dist * direction.z() < 0)
return std::nullopt;
102 const double direction_scale = dist / direction.z();
103 return start + direction * direction_scale;
107 std::optional<TVector3> cylinder_line_intersection(TVector3
s, TVector3
v,
double radius)
114 const double b = 2 * (v.x() * s.x() + v.y() * s.y());
117 const double rootterm =
square(b) - 4 * a *
c;
124 if (rootterm < 0 || a == 0)
return std::nullopt;
127 const double sqrtterm = std::sqrt(rootterm);
128 const double t1 = (-b + sqrtterm) / (2 * a);
129 const double t2 = (-b - sqrtterm) / (2 * a);
135 const double& min_t = (t2 < t1 && t2 > 0) ? t2 : t1;
136 return s + v * min_t;
140 std::optional<TVector3> field_cage_intersection(TVector3 start, TVector3 direction)
142 const auto ofc_strike = cylinder_line_intersection(start, direction, end_CM);
143 const auto ifc_strike = cylinder_line_intersection(start, direction, begin_CM);
146 if (!ifc_strike)
return ofc_strike;
147 if (!ofc_strike)
return ifc_strike;
150 const auto ifc_dist = (ifc_strike->Z() - start.Z()) / direction.Z();
151 const auto ofc_dist = (ofc_strike->Z() - start.Z()) / direction.Z();
154 return (ofc_dist > 0) ? ofc_strike : std::nullopt;
155 else if (ofc_dist < 0)
158 return (ifc_dist < ofc_dist) ? ifc_strike : ofc_strike;
162 inline std::ostream&
operator<<(std::ostream&
out,
const TVector3& vector)
164 out <<
"( " << vector.x() <<
", " << vector.y() <<
", " << vector.z() <<
")";
184 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
187 std::cout <<
"Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
193 std::cout <<
PHWHERE <<
"DST Node missing, aborting." << std::endl;
203 auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode,
hitnodename.c_str());
206 std::cout <<
Name() <<
" Could not locate G4HIT node " <<
hitnodename << std::endl;
220 std::cout <<
PHWHERE <<
"DST Node missing, aborting." << std::endl;
245 std::cout<<
"PHG4TpcDirectLaser::InitRun - m_steppingpattern: " <<
m_steppingpattern << std::endl;
246 std::cout<<
"PHG4TpcDirectLaser::InitRun - nTotalSteps: " <<
nTotalSteps << std::endl;
251 std::cout <<
"PHG4TpcDirectLaser::InitRun - phi steps: " <<
nPhiSteps <<
" min: " <<
minPhi <<
" max: " <<
maxPhi << std::endl;
252 std::cout <<
"PHG4TpcDirectLaser::InitRun - theta steps: " <<
nThetaSteps <<
" min: " <<
minTheta <<
" max: " <<
maxTheta << std::endl;
253 std::cout <<
"PHG4TpcDirectLaser::InitRun - nTotalSteps: " <<
nTotalSteps << std::endl;
255 std::cout <<
"PHG4TpcDirectLaser::InitRun - electrons_per_cm: " <<
electrons_per_cm << std::endl;
256 std::cout <<
"PHG4TpcDirectLaser::InitRun - electrons_per_gev " <<
electrons_per_gev << std::endl;
260 std::string LASER_ANGLES_ROOTFILE =
std::string(getenv(
"CALIBRATIONROOT")) +
"/TPC/DirectLaser/theta_phi_laser.root";
261 TFile* infile1 = TFile::Open(LASER_ANGLES_ROOTFILE.c_str());
263 pattern = (TNtuple*) infile1->Get(
"angles");
264 pattern->SetBranchAddress(
"#theta",&
theta_p);
265 pattern->SetBranchAddress(
"#phi",&
phi_p);
274 m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,
"G4TruthInfo");
315 static constexpr
double Ne_dEdx = 1.56;
316 static constexpr
double CF4_dEdx = 7.00;
319 static constexpr
double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 *
CF4_NTotal;
320 static constexpr
double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 *
CF4_dEdx;
321 static constexpr
double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
333 if (n < 0 || max < min)
335 std::cout <<
PHWHERE <<
" - invalid" << std::endl;
347 if (n < 0 || max < min)
349 std::cout <<
PHWHERE <<
" - invalid" << std::endl;
363 if (n < 0 || n > 13802)
365 std::cout <<
PHWHERE <<
" - invalid" << std::endl;
381 const TVector3 position_base(60 * cm, 0., halflength_tpc);
384 for (
int i = 0;
i < 8; ++
i)
398 laser.
m_phi = M_PI / 2 *
i + (15 * M_PI/180);
404 laser.
m_phi = M_PI / 2 *
i - (15 * M_PI/180);
441 std::cout <<
"PHG4TpcDirectLaser::AimToThetaPhi - theta: " << theta <<
" phi: " << phi << std::endl;
459 std::cout <<
"PHG4TpcDirectLaser::AimToPatternStep - step: " << n <<
"/" <<
nTotalSteps << std::endl;
488 std::cout <<
"PHG4TpcDirectLaser::AimToPatternStep_File - step: " << n <<
"/" <<
nTotalSteps << std::endl;
497 std::cout <<
"From file, current entry = " << n <<
" Theta: " <<
theta_p <<
" Phi: " <<
phi_p << std::endl;
518 std::cout <<
PHWHERE <<
"invalid g4hit container. aborting" << std::endl;
527 TVector3 dir(0, 0, direction);
530 dir.RotateY(theta * direction);
533 else dir.RotateZ(-phi);
536 dir.RotateZ(laser.
m_phi);
541 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack - position: " <<
pos <<
" direction: " << dir << std::endl;
545 static constexpr
double total_momentum = 1;
555 const auto vertex =
new PHG4VtxPoint_t(
pos.x(),
pos.y(),
pos.z(), 0, vtxid);
562 auto particle =
new PHG4Particle_t();
567 particle->set_px(total_momentum * dir.x());
568 particle->set_py(total_momentum * dir.y());
569 particle->set_pz(total_momentum * dir.z());
583 track.
set_px(total_momentum * dir.x());
584 track.
set_py(total_momentum * dir.y());
585 track.
set_pz(total_momentum * dir.z());
592 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack - position: " <<
pos <<
" direction: " << dir << std::endl;
602 const auto plane_strike = (
pos.z() * dir.z() > 0) ? endcap_intersection(
pos, dir) : central_membrane_intersection(
pos, dir);
605 const auto fc_strike = field_cage_intersection(
pos, dir);
608 if (!(plane_strike || fc_strike))
return;
612 const TVector3& strike = (fc_strike && (!plane_strike || fc_strike->z() / dir.z() < plane_strike->z() / dir.z())) ? *fc_strike : *plane_strike;
616 double fullLength = delta.Mag();
617 int nHitSteps = fullLength / maxHitLength + 1;
619 TVector3 start =
pos;
620 TVector3 end =
start;
621 TVector3
step = dir * (maxHitLength / (dir.Mag()));
622 double stepLength = 0;
626 std::cout <<
"PHG4TpcDirectLaser::AppendLaserTrack -"
627 <<
" fullLength: " << fullLength
628 <<
" nHitSteps: " << nHitSteps
632 for (
int i = 0;
i < nHitSteps;
i++)
635 if (
i + 1 == nHitSteps)
640 stepLength = delta.Mag();
646 stepLength = step.Mag();
650 auto hit =
new PHG4Hit_t;
651 hit->set_trkid(trackid);
655 hit->set_x(0, start.X() /
cm);
656 hit->set_y(0, start.Y() /
cm);
657 hit->set_z(0, start.Z() /
cm);
658 hit->set_t(0, (start -
pos).Mag() / speed_of_light);
660 hit->set_x(1, end.X() /
cm);
661 hit->set_y(1, end.Y() /
cm);
662 hit->set_z(1, end.Z() /
cm);
663 hit->set_t(1, (end -
pos).Mag() / speed_of_light);
666 hit->set_px(0, dir.X());
667 hit->set_py(0, dir.Y());
668 hit->set_pz(0, dir.Z());
670 hit->set_px(1, dir.X());
671 hit->set_py(1, dir.Y());
672 hit->set_pz(1, dir.Z());
676 hit->set_eion(totalE);
677 hit->set_edep(totalE);