ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcDirectLaser.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcDirectLaser.cc
1 #include "PHG4TpcDirectLaser.h"
2 
3 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
4 
6 #include <g4main/PHG4HitDefs.h> // for get_volume_id
7 #include <g4main/PHG4Hitv1.h>
10 #include <g4main/PHG4VtxPointv1.h>
11 
13 #include <fun4all/SubsysReco.h> // for SubsysReco
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h>
17 #include <phool/PHNode.h>
18 #include <phool/PHNodeIterator.h>
19 #include <phool/PHObject.h> // for PHObject
20 #include <phool/getClass.h>
21 #include <phool/phool.h> // for PHWHERE
22 
26 
27 #include <TVector3.h> // for TVector3, operator*
28 
29 #include <gsl/gsl_const_mksa.h> // for the speed of light
30 
31 #include <cassert>
32 #include <iostream> // for operator<<, basic_os...
33 #include <optional>
34 
35 namespace
36 {
37  using PHG4Particle_t = PHG4Particlev3;
38  using PHG4VtxPoint_t = PHG4VtxPointv1;
39  using PHG4Hit_t = PHG4Hitv1;
40 
41  // utility
42  template <class T>
43  inline constexpr T square(const T& x)
44  {
45  return x * x;
46  }
47 
48  // unique detector id for all direct lasers
49  static const int detId = PHG4HitDefs::get_volume_id("PHG4TpcDirectLaser");
50 
52 
53  static constexpr double cm = 1.0;
55 
57  static constexpr double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1e-7;
58 
60  static constexpr double maxHitLength = 1. * cm;
61 
63  static constexpr double halflength_tpc = 105.5 * cm;
64 
65  // inner and outer radii of field cages/TPC
66  static constexpr double begin_CM = 20. * cm;
67  static constexpr double end_CM = 78. * cm;
68 
69  //half the thickness of the CM;
70  static constexpr double halfwidth_CM = 0.5 * cm;
71 
72  //_____________________________________________________________
73  std::optional<TVector3> central_membrane_intersection(TVector3 start, TVector3 direction)
74  {
75  const double end = start.z() > 0 ? halfwidth_CM : -halfwidth_CM;
76  const double dist = end - start.z();
77 
78  // if line is vertical, it will never intercept the endcap
79  if (direction.z() == 0) return std::nullopt;
80 
81  // check that distance and direction have the same sign
82  if (dist * direction.z() < 0) return std::nullopt;
83 
84  const double direction_scale = dist / direction.z();
85  return start + direction * direction_scale;
86  }
87 
88  //_____________________________________________________________
89  std::optional<TVector3> endcap_intersection(TVector3 start, TVector3 direction)
90  {
91  const double end = start.z() > 0 ? halflength_tpc : -halflength_tpc;
92  const double dist = end - start.z();
93 
94  // if line is vertical, it will never intercept the endcap
95  if (direction.z() == 0) return std::nullopt;
96 
97  // check that distance and direction have the same sign
98  if (dist * direction.z() < 0) return std::nullopt;
99 
100  const double direction_scale = dist / direction.z();
101  return start + direction * direction_scale;
102  }
103 
104  //_____________________________________________________________
105  std::optional<TVector3> cylinder_line_intersection(TVector3 s, TVector3 v, double radius)
106  {
107  const double R2 = square(radius);
108 
109  //Generalized Parameters for collision with cylinder of radius R:
110  //from quadratic formula solutions of when a vector intersects a circle:
111  const double a = square(v.x()) + square(v.y());
112  const double b = 2 * (v.x() * s.x() + v.y() * s.y());
113  const double c = square(s.x()) + square(s.y()) - R2;
114 
115  const double rootterm = square(b) - 4 * a * c;
116 
117  /*
118  * if a==0 then we are parallel and will have no solutions.
119  * if the rootterm is negative, we will have no real roots,
120  * we are outside the cylinder and pointing skew to the cylinder such that we never cross.
121  */
122  if (rootterm < 0 || a == 0) return std::nullopt;
123 
124  //Find the (up to) two points where we collide with the cylinder:
125  const double sqrtterm = std::sqrt(rootterm);
126  const double t1 = (-b + sqrtterm) / (2 * a);
127  const double t2 = (-b - sqrtterm) / (2 * a);
128 
129  /*
130  * if either of the t's are nonzero, we have a collision
131  * the collision closest to the start (hence with the smallest t that is greater than zero) is the one that happens.
132  */
133  const double& min_t = (t2 < t1 && t2 > 0) ? t2 : t1;
134  return s + v * min_t;
135  }
136 
137  //_____________________________________________________________
138  std::optional<TVector3> field_cage_intersection(TVector3 start, TVector3 direction)
139  {
140  const auto ofc_strike = cylinder_line_intersection(start, direction, end_CM);
141  const auto ifc_strike = cylinder_line_intersection(start, direction, begin_CM);
142 
143  // if either of the two intersection is invalid, return the other
144  if (!ifc_strike) return ofc_strike;
145  if (!ofc_strike) return ifc_strike;
146 
147  // both intersection are valid, calculate signed distance to start z
148  const auto ifc_dist = (ifc_strike->Z() - start.Z()) / direction.Z();
149  const auto ofc_dist = (ofc_strike->Z() - start.Z()) / direction.Z();
150 
151  if (ifc_dist < 0)
152  return (ofc_dist > 0) ? ofc_strike : std::nullopt;
153  else if (ofc_dist < 0)
154  return ifc_strike;
155  else
156  return (ifc_dist < ofc_dist) ? ifc_strike : ofc_strike;
157  }
158 
160  inline std::ostream& operator<<(std::ostream& out, const TVector3& vector)
161  {
162  out << "( " << vector.x() << ", " << vector.y() << ", " << vector.z() << ")";
163  return out;
164  }
165 
166 } // namespace
167 
168 //_____________________________________________________________
170  : SubsysReco(name)
171  , PHParameterInterface(name)
172 {
174 }
175 
176 //_____________________________________________________________
178 {
179  // g4 truth info
180  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
181  if (!m_g4truthinfo)
182  {
183  std::cout << "Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
184 
185  PHNodeIterator iter(topNode);
186  auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
187  if (!dstNode)
188  {
189  std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
191  }
192 
194  dstNode->addNode(new PHIODataNode<PHObject>(m_g4truthinfo, "G4TruthInfo", "PHObject"));
195  }
196 
197  // load and check G4Hit node
198  hitnodename = "G4HIT_" + detector;
199  auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
200  if (!g4hit)
201  {
202  std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
204  }
205 
206  // find or create track map
207  /* it is used to store laser parameters on a per event basis */
208  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
209  if (!m_track_map)
210  {
211  // find DST node and check
212  PHNodeIterator iter(topNode);
213  auto dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
214  if (!dstNode)
215  {
216  std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
218  }
219 
220  // find or create SVTX node
221  iter = PHNodeIterator(dstNode);
222  auto node = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "SVTX"));
223  if (!node) dstNode->addNode(node = new PHCompositeNode("SVTX"));
224 
225  // add track node
227  node->addNode(new PHIODataNode<PHObject>(m_track_map, m_track_map_name, "PHObject"));
228  }
229 
230  // setup parameters
232  electrons_per_cm = get_int_param("electrons_per_cm");
233  electrons_per_gev = get_double_param("electrons_per_gev");
234 
235  // setup lasers
236  SetupLasers();
237 
238  // print configuration
239  std::cout << "PHG4TpcDirectLaser::InitRun - m_autoAdvanceDirectLaser: " << m_autoAdvanceDirectLaser << 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;
243 
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;
246 
248 }
249 
250 //_____________________________________________________________
252 {
253  // g4 input event
254  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
255  assert(m_g4truthinfo);
256 
257  // load g4hit container
258  m_g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
259  assert(m_g4hitcontainer);
260 
261  // load track map
262  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
263  assert(m_track_map);
264 
266  {
268  }
269  else
270  {
271  // use arbitrary direction
272  AimToThetaPhi(M_PI / 180. * 10., M_PI / 180. * 90);
273  }
274 
276 }
277 
278 //_____________________________________________________________
280 {
281  // same gas parameters as in PHG4TpcElectronDrift::SetDefaultParameters
282 
283  // Data on gasses @20 C and 760 Torr from the following source:
284  // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
285  // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
286  // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
287  static constexpr double Ne_dEdx = 1.56; // keV/cm
288  static constexpr double CF4_dEdx = 7.00; // keV/cm
289  static constexpr double Ne_NTotal = 43; // Number/cm
290  static constexpr double CF4_NTotal = 100; // Number/cm
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;
294 
295  // number of electrons per deposited GeV in TPC gas
296  set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1e6);
297 
298  // number of electrons deposited by laser per cm
299  set_default_int_param("electrons_per_cm", 72);
300 }
301 
302 //_____________________________________________________________
303 void PHG4TpcDirectLaser::SetPhiStepping(int n, double min, double max)
304 {
305  if (n < 0 || max < min)
306  {
307  std::cout << PHWHERE << " - invalid" << std::endl;
308  return;
309  }
310  nPhiSteps = n;
311  minPhi = min;
312  maxPhi = max;
314  return;
315 }
316 //_____________________________________________________________
318 {
319  if (n < 0 || max < min)
320  {
321  std::cout << PHWHERE << " - invalid" << std::endl;
322  return;
323  }
324  nThetaSteps = n;
325  minTheta = min;
326  maxTheta = max;
328 
329  return;
330 }
331 
332 //_____________________________________________________________
334 {
335  // clear previous lasers
336  m_lasers.clear();
337 
338  // position of first laser at positive z
339  const TVector3 position_base(60 * cm, 0., halflength_tpc);
340 
341  // add lasers
342  for (int i = 0; i < 8; ++i)
343  {
344  Laser laser;
345 
346  // set laser direction
347  /*
348  * first four lasers are on positive z readout plane, and shoot towards negative z
349  * next four lasers are on negative z readout plane and shoot towards positive z
350  */
351  laser.m_position = position_base;
352  if (i < 4)
353  {
354  laser.m_position.SetZ(position_base.z());
355  laser.m_direction = -1;
356  }
357  else
358  {
359  laser.m_position.SetZ(-position_base.z());
360  laser.m_direction = 1;
361  }
362 
363  // rotate around z
364  laser.m_phi = M_PI / 2 * i;
365  laser.m_position.RotateZ(laser.m_phi);
366 
367  // append
368  m_lasers.push_back(laser);
369  }
370 }
371 
372 //_____________________________________________________________
374 {
375  if (nTotalSteps >= 1)
376  {
379  }
380 }
381 
382 //_____________________________________________________________
384 {
385  if (Verbosity())
386  {
387  std::cout << "PHG4TpcDirectLaser::AimToThetaPhi - theta: " << theta << " phi: " << phi << std::endl;
388  }
389 
390  // all lasers
391  for (const auto& laser : m_lasers)
392  {
393  AppendLaserTrack(theta, phi, laser);
394  }
395 }
396 
397 //_____________________________________________________________
399 {
400  //trim against overflows
401  n = n % nTotalSteps;
402 
403  if (Verbosity())
404  {
405  std::cout << "PHG4TpcDirectLaser::AimToPatternStep - step: " << n << "/" << nTotalSteps << std::endl;
406  }
407 
408  // store as current pattern
410 
411  // calculate theta
412  const int thetaStep = n / nPhiSteps;
413  const double theta = minTheta + thetaStep * (maxTheta - minTheta) / nThetaSteps;
414 
415  // calculate phi
416  const int phiStep = n % nPhiSteps;
417  const double phi = minPhi + phiStep * (maxPhi - minPhi) / nPhiSteps;
418 
419  // generate laser tracks
420  AimToThetaPhi(theta, phi);
421 
422  return;
423 }
424 
425 //_____________________________________________________________
427 {
428  if (!m_g4hitcontainer)
429  {
430  std::cout << PHWHERE << "invalid g4hit container. aborting" << std::endl;
431  return;
432  }
433 
434  // store laser position
435  const auto& pos = laser.m_position;
436 
437  // define track direction
438  const auto& direction = laser.m_direction;
439  TVector3 dir(0, 0, direction);
440 
441  //adjust direction
442  dir.RotateY(theta * direction);
443  dir.RotateZ(phi);
444 
445  // also rotate by laser azimuth
446  dir.RotateZ(laser.m_phi);
447 
448  // print
449  if (Verbosity())
450  {
451  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
452  }
453 
454  // dummy momentum
455  static constexpr double total_momentum = 1;
456 
457  // mc track id
458  int trackid = -1;
459 
460  // create truth vertex and particle
461  if (m_g4truthinfo)
462  {
463  // add vertex
464  const auto vtxid = m_g4truthinfo->maxvtxindex() + 1;
465  const auto vertex = new PHG4VtxPoint_t(pos.x(), pos.y(), pos.z(), 0, vtxid);
466  m_g4truthinfo->AddVertex(vtxid, vertex);
467 
468  // increment track id
469  trackid = m_g4truthinfo->maxtrkindex() + 1;
470 
471  // create new g4particle
472  auto particle = new PHG4Particle_t();
473  particle->set_track_id(trackid);
474  particle->set_vtx_id(vtxid);
475  particle->set_parent_id(0);
476  particle->set_primary_id(trackid);
477  particle->set_px(total_momentum * dir.x());
478  particle->set_py(total_momentum * dir.y());
479  particle->set_pz(total_momentum * dir.z());
480 
482  }
483 
484  // store in SvtxTrack map
485  if (m_track_map)
486  {
488  track.set_x(pos.x());
489  track.set_y(pos.y());
490  track.set_z(pos.z());
491 
492  // total momentum is irrelevant. What matters is the direction
493  track.set_px(total_momentum * dir.x());
494  track.set_py(total_momentum * dir.y());
495  track.set_pz(total_momentum * dir.z());
496 
497  // insert in map
498  m_track_map->insert(&track);
499 
500  if (Verbosity())
501  {
502  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
503  }
504  }
505 
506  // find collision point
507  /*
508  * intersection to either central membrane or endcaps
509  * if the position along beam and laser direction have the same sign, it will intercept the endcap
510  * otherwise will intercept the central membrane
511  */
512  const auto plane_strike = (pos.z() * dir.z() > 0) ? endcap_intersection(pos, dir) : central_membrane_intersection(pos, dir);
513 
514  // field cage intersection
515  const auto fc_strike = field_cage_intersection(pos, dir);
516 
517  // if none of the strikes is valid, there is no valid information found.
518  if (!(plane_strike || fc_strike)) return;
519 
520  // decide relevant end of laser
521  /* chose field cage intersection if valid, and if either plane intersection is invalid or happens on a larger z along the laser direction) */
522  const TVector3& strike = (fc_strike && (!plane_strike || fc_strike->z() / dir.z() < plane_strike->z() / dir.z())) ? *fc_strike : *plane_strike;
523 
524  //find length
525  TVector3 delta = (strike - pos);
526  double fullLength = delta.Mag();
527  int nHitSteps = fullLength / maxHitLength + 1;
528 
529  TVector3 start = pos;
530  TVector3 end = start;
531  TVector3 step = dir * (maxHitLength / (dir.Mag()));
532  double stepLength = 0;
533 
534  if (Verbosity())
535  {
536  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack -"
537  << " fullLength: " << fullLength
538  << " nHitSteps: " << nHitSteps
539  << std::endl;
540  }
541 
542  for (int i = 0; i < nHitSteps; i++)
543  {
544  start = end; //new starting point is the previous ending point.
545  if (i + 1 == nHitSteps)
546  {
547  //last step is the remainder size
548  end = strike;
549  delta = start - end;
550  stepLength = delta.Mag();
551  }
552  else
553  {
554  //all other steps are uniform length
555  end = start + step;
556  stepLength = step.Mag();
557  }
558 
559  //from phg4tpcsteppingaction.cc
560  auto hit = new PHG4Hit_t;
561  hit->set_trkid(trackid);
562  hit->set_layer(99);
563 
564  //here we set the entrance values in cm
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);
569 
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);
574 
575  // momentum
576  hit->set_px(0, dir.X()); // GeV
577  hit->set_py(0, dir.Y());
578  hit->set_pz(0, dir.Z());
579 
580  hit->set_px(1, dir.X());
581  hit->set_py(1, dir.Y());
582  hit->set_pz(1, dir.Z());
583 
584  const double totalE = electrons_per_cm * stepLength / electrons_per_gev;
585 
586  hit->set_eion(totalE);
587  hit->set_edep(totalE);
588  m_g4hitcontainer->AddHit(detId, hit);
589  }
590 
591  return;
592 }