Analysis Software
Documentation for sPHENIX simulation 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 #include <TNtuple.h>
29 #include <TFile.h>
30 
31 #include <gsl/gsl_const_mksa.h> // for the speed of light
32 
33 #include <cassert>
34 #include <iostream> // for operator<<, basic_os...
35 #include <optional>
36 
37 namespace
38 {
39  using PHG4Particle_t = PHG4Particlev3;
40  using PHG4VtxPoint_t = PHG4VtxPointv1;
41  using PHG4Hit_t = PHG4Hitv1;
42 
43  // utility
44  template <class T>
45  inline constexpr T square(const T& x)
46  {
47  return x * x;
48  }
49 
50  // unique detector id for all direct lasers
51  static const int detId = PHG4HitDefs::get_volume_id("PHG4TpcDirectLaser");
52 
54 
55  static constexpr double cm = 1.0;
57 
59  static constexpr double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1e-7;
60 
62  static constexpr double maxHitLength = 1. * cm;
63 
65  static constexpr double halflength_tpc = 105.5 * cm;
66 
67  // inner and outer radii of field cages/TPC
68  static constexpr double begin_CM = 20. * cm;
69  static constexpr double end_CM = 78. * cm;
70 
71  //half the thickness of the CM;
72  static constexpr double halfwidth_CM = 0.5 * cm;
73 
74  //_____________________________________________________________
75  std::optional<TVector3> central_membrane_intersection(TVector3 start, TVector3 direction)
76  {
77  const double end = start.z() > 0 ? halfwidth_CM : -halfwidth_CM;
78  const double dist = end - start.z();
79 
80  // if line is vertical, it will never intercept the endcap
81  if (direction.z() == 0) return std::nullopt;
82 
83  // check that distance and direction have the same sign
84  if (dist * direction.z() < 0) return std::nullopt;
85 
86  const double direction_scale = dist / direction.z();
87  return start + direction * direction_scale;
88  }
89 
90  //_____________________________________________________________
91  std::optional<TVector3> endcap_intersection(TVector3 start, TVector3 direction)
92  {
93  const double end = start.z() > 0 ? halflength_tpc : -halflength_tpc;
94  const double dist = end - start.z();
95 
96  // if line is vertical, it will never intercept the endcap
97  if (direction.z() == 0) return std::nullopt;
98 
99  // check that distance and direction have the same sign
100  if (dist * direction.z() < 0) return std::nullopt;
101 
102  const double direction_scale = dist / direction.z();
103  return start + direction * direction_scale;
104  }
105 
106  //_____________________________________________________________
107  std::optional<TVector3> cylinder_line_intersection(TVector3 s, TVector3 v, double radius)
108  {
109  const double R2 = square(radius);
110 
111  //Generalized Parameters for collision with cylinder of radius R:
112  //from quadratic formula solutions of when a vector intersects a circle:
113  const double a = square(v.x()) + square(v.y());
114  const double b = 2 * (v.x() * s.x() + v.y() * s.y());
115  const double c = square(s.x()) + square(s.y()) - R2;
116 
117  const double rootterm = square(b) - 4 * a * c;
118 
119  /*
120  * if a==0 then we are parallel and will have no solutions.
121  * if the rootterm is negative, we will have no real roots,
122  * we are outside the cylinder and pointing skew to the cylinder such that we never cross.
123  */
124  if (rootterm < 0 || a == 0) return std::nullopt;
125 
126  //Find the (up to) two points where we collide with the cylinder:
127  const double sqrtterm = std::sqrt(rootterm);
128  const double t1 = (-b + sqrtterm) / (2 * a);
129  const double t2 = (-b - sqrtterm) / (2 * a);
130 
131  /*
132  * if either of the t's are nonzero, we have a collision
133  * the collision closest to the start (hence with the smallest t that is greater than zero) is the one that happens.
134  */
135  const double& min_t = (t2 < t1 && t2 > 0) ? t2 : t1;
136  return s + v * min_t;
137  }
138 
139  //_____________________________________________________________
140  std::optional<TVector3> field_cage_intersection(TVector3 start, TVector3 direction)
141  {
142  const auto ofc_strike = cylinder_line_intersection(start, direction, end_CM);
143  const auto ifc_strike = cylinder_line_intersection(start, direction, begin_CM);
144 
145  // if either of the two intersection is invalid, return the other
146  if (!ifc_strike) return ofc_strike;
147  if (!ofc_strike) return ifc_strike;
148 
149  // both intersection are valid, calculate signed distance to start z
150  const auto ifc_dist = (ifc_strike->Z() - start.Z()) / direction.Z();
151  const auto ofc_dist = (ofc_strike->Z() - start.Z()) / direction.Z();
152 
153  if (ifc_dist < 0)
154  return (ofc_dist > 0) ? ofc_strike : std::nullopt;
155  else if (ofc_dist < 0)
156  return ifc_strike;
157  else
158  return (ifc_dist < ofc_dist) ? ifc_strike : ofc_strike;
159  }
160 
162  inline std::ostream& operator<<(std::ostream& out, const TVector3& vector)
163  {
164  out << "( " << vector.x() << ", " << vector.y() << ", " << vector.z() << ")";
165  return out;
166  }
167 
168 } // namespace
169 
170 //_____________________________________________________________
172  : SubsysReco(name)
173  , PHParameterInterface(name)
174 {
176  theta_p = 0;
177  phi_p = 0;
178 }
179 
180 //_____________________________________________________________
182 {
183  // g4 truth info
184  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
185  if (!m_g4truthinfo)
186  {
187  std::cout << "Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
188 
189  PHNodeIterator iter(topNode);
190  auto dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
191  if (!dstNode)
192  {
193  std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
195  }
196 
198  dstNode->addNode(new PHIODataNode<PHObject>(m_g4truthinfo, "G4TruthInfo", "PHObject"));
199  }
200 
201  // load and check G4Hit node
202  hitnodename = "G4HIT_" + detector;
203  auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
204  if (!g4hit)
205  {
206  std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
208  }
209 
210  // find or create track map
211  /* it is used to store laser parameters on a per event basis */
212  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
213  if (!m_track_map)
214  {
215  // find DST node and check
216  PHNodeIterator iter(topNode);
217  auto dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
218  if (!dstNode)
219  {
220  std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
222  }
223 
224  // find or create SVTX node
225  iter = PHNodeIterator(dstNode);
226  auto node = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "SVTX"));
227  if (!node) dstNode->addNode(node = new PHCompositeNode("SVTX"));
228 
229  // add track node
231  node->addNode(new PHIODataNode<PHObject>(m_track_map, m_track_map_name, "PHObject"));
232  }
233 
234  // setup parameters
236  electrons_per_cm = get_int_param("electrons_per_cm");
237  electrons_per_gev = get_double_param("electrons_per_gev");
238 
239  // setup lasers
240  SetupLasers();
241 
242  // print configuration
243  if(m_steppingpattern == true)
244  {
245  std::cout<< "PHG4TpcDirectLaser::InitRun - m_steppingpattern: " << m_steppingpattern << std::endl;
246  std::cout<< "PHG4TpcDirectLaser::InitRun - nTotalSteps: " << nTotalSteps << std::endl;
247  }
248  else
249  {
250  std::cout << "PHG4TpcDirectLaser::InitRun - m_autoAdvanceDirectLaser: " << m_autoAdvanceDirectLaser << 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;
254  }
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;
257 
258  //TFile * infile1 = TFile::Open("theta_phi_laser.root");
259 
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());
262 
263  pattern = (TNtuple*) infile1->Get("angles");
264  pattern->SetBranchAddress("#theta",&theta_p);
265  pattern->SetBranchAddress("#phi",&phi_p);
266 
268 }
269 
270 //_____________________________________________________________
272 {
273  // g4 input event
274  m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
276 
277  // load g4hit container
278  m_g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
280 
281  // load track map
282  m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
284 
286  {
288  }
289  //_________________________________________________
290 
291  else if (m_steppingpattern)
292  {
294  }
295 
296  //_________________________________________________
297  else
298  {
299  // use arbitrary direction
301  }
302 
304 }
305 
306 //_____________________________________________________________
308 {
309  // same gas parameters as in PHG4TpcElectronDrift::SetDefaultParameters
310 
311  // Data on gasses @20 C and 760 Torr from the following source:
312  // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
313  // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
314  // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
315  static constexpr double Ne_dEdx = 1.56; // keV/cm
316  static constexpr double CF4_dEdx = 7.00; // keV/cm
317  static constexpr double Ne_NTotal = 43; // Number/cm
318  static constexpr double CF4_NTotal = 100; // Number/cm
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;
322 
323  // number of electrons per deposited GeV in TPC gas
324  set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1e6);
325 
326  // number of electrons deposited by laser per cm
327  set_default_int_param("electrons_per_cm", 72);
328 }
329 
330 //_____________________________________________________________
331 void PHG4TpcDirectLaser::SetPhiStepping(int n, double min, double max)
332 {
333  if (n < 0 || max < min)
334  {
335  std::cout << PHWHERE << " - invalid" << std::endl;
336  return;
337  }
338  nPhiSteps = n;
339  minPhi = min;
340  maxPhi = max;
342  return;
343 }
344 //_____________________________________________________________
345 void PHG4TpcDirectLaser::SetThetaStepping(int n, double min, double max)
346 {
347  if (n < 0 || max < min)
348  {
349  std::cout << PHWHERE << " - invalid" << std::endl;
350  return;
351  }
352  nThetaSteps = n;
353  minTheta = min;
354  maxTheta = max;
356 
357  return;
358 }
359 
360 //_____________________________________________________________
362 {
363  if (n < 0 || n > 13802) //13802 = hard coded number of tuple entries
364  {
365  std::cout << PHWHERE << " - invalid" << std::endl;
366  return;
367  }
368  nTotalSteps = n;
369 
370  return;
371 }
372 
373 //_____________________________________________________________
374 
376 {
377  // clear previous lasers
378  m_lasers.clear();
379 
380  // position of first laser at positive z
381  const TVector3 position_base(60 * cm, 0., halflength_tpc);
382 
383  // add lasers
384  for (int i = 0; i < 8; ++i)
385  {
386  Laser laser;
387 
388  // set laser direction
389  /*
390  * first four lasers are on positive z readout plane, and shoot towards negative z
391  * next four lasers are on negative z readout plane and shoot towards positive z
392  */
393  laser.m_position = position_base;
394  if (i < 4)
395  {
396  laser.m_position.SetZ(position_base.z());
397  laser.m_direction = -1;
398  laser.m_phi = M_PI / 2 * i + (15 * M_PI/180); //additional offset of 15 deg.
399  }
400  else
401  {
402  laser.m_position.SetZ(-position_base.z());
403  laser.m_direction = 1;
404  laser.m_phi = M_PI / 2 * i - (15 * M_PI/180); //additional offset of 15 deg.
405  }
406 
407  // rotate around z
408  laser.m_position.RotateZ(laser.m_phi);
409 
410  // append
411  m_lasers.push_back(laser); //All lasers
412  // if(i==0) m_lasers.push_back(laser);//Only laser 1
413  // if(i==3) m_lasers.push_back(laser);// Laser 4
414  // if(i<4) m_lasers.push_back(laser);//Lasers 1, 2, 3, 4
415  }
416 }
417 
418 //_____________________________________________________________
420 {
421  if (nTotalSteps >= 1)
422  {
423  if (m_steppingpattern)
424  {
427  }
428  else
429  {
432  }
433  }
434 }
435 
436 //_____________________________________________________________
438 {
439  if (Verbosity())
440  {
441  std::cout << "PHG4TpcDirectLaser::AimToThetaPhi - theta: " << theta << " phi: " << phi << std::endl;
442  }
443 
444  // all lasers
445  for (const auto& laser : m_lasers)
446  {
447  AppendLaserTrack(theta, phi, laser);
448  }
449 }
450 
451 //_____________________________________________________________
453 {
454  //trim against overflows
455  n = n % nTotalSteps;
456 
457  if (Verbosity())
458  {
459  std::cout << "PHG4TpcDirectLaser::AimToPatternStep - step: " << n << "/" << nTotalSteps << std::endl;
460  }
461 
462  // store as current pattern
464 
465  // calculate theta
466  const int thetaStep = n / nPhiSteps;
467  const double theta = minTheta + thetaStep * (maxTheta - minTheta) / nThetaSteps;
468 
469  // calculate phi
470  const int phiStep = n % nPhiSteps;
471  const double phi = minPhi + phiStep * (maxPhi - minPhi) / nPhiSteps;
472 
473  // generate laser tracks
474  AimToThetaPhi(theta, phi);
475 
476  return;
477 }
478 
479 //_____________________________________________________________
480 
482 {
483  //trim against overflows
484  n = n % nTotalSteps;
485 
486  if (Verbosity())
487  {
488  std::cout << "PHG4TpcDirectLaser::AimToPatternStep_File - step: " << n << "/" << nTotalSteps << std::endl;
489  }
490 
491  // store as current pattern
493 
494  pattern->GetEntry(n);
495 
496  // calculate theta
497  std::cout << "From file, current entry = " << n << " Theta: " << theta_p <<" Phi: " << phi_p << std::endl;
498 
499  const double theta = theta_p*M_PI/180.;
500 
501  // calculate phi
502  const double phi = phi_p*M_PI/180.;
503 
504 
505  // generate laser tracks
506  AimToThetaPhi(theta, phi);
507 
508  return;
509 }
510 
511 //_____________________________________________________________
512 
513 
515 {
516  if (!m_g4hitcontainer)
517  {
518  std::cout << PHWHERE << "invalid g4hit container. aborting" << std::endl;
519  return;
520  }
521 
522  // store laser position
523  const auto& pos = laser.m_position;
524 
525  // define track direction
526  const auto& direction = laser.m_direction;
527  TVector3 dir(0, 0, direction);
528 
529  //adjust direction
530  dir.RotateY(theta * direction);
531 
532  if(laser.m_direction == -1) dir.RotateZ(phi); //if +z facing -z
533  else dir.RotateZ(-phi); //if -z facting +z
534 
535  // also rotate by laser azimuth
536  dir.RotateZ(laser.m_phi);
537 
538  // print
539  if (Verbosity())
540  {
541  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
542  }
543 
544  // dummy momentum
545  static constexpr double total_momentum = 1;
546 
547  // mc track id
548  int trackid = -1;
549 
550  // create truth vertex and particle
551  if (m_g4truthinfo)
552  {
553  // add vertex
554  const auto vtxid = m_g4truthinfo->maxvtxindex() + 1;
555  const auto vertex = new PHG4VtxPoint_t(pos.x(), pos.y(), pos.z(), 0, vtxid);
556  m_g4truthinfo->AddVertex(vtxid, vertex);
557 
558  // increment track id
559  trackid = m_g4truthinfo->maxtrkindex() + 1;
560 
561  // create new g4particle
562  auto particle = new PHG4Particle_t();
563  particle->set_track_id(trackid);
564  particle->set_vtx_id(vtxid);
565  particle->set_parent_id(0);
566  particle->set_primary_id(trackid);
567  particle->set_px(total_momentum * dir.x());
568  particle->set_py(total_momentum * dir.y());
569  particle->set_pz(total_momentum * dir.z());
570 
572  }
573 
574  // store in SvtxTrack map
575  if (m_track_map)
576  {
577  SvtxTrack_v2 track;
578  track.set_x(pos.x());
579  track.set_y(pos.y());
580  track.set_z(pos.z());
581 
582  // total momentum is irrelevant. What matters is the direction
583  track.set_px(total_momentum * dir.x());
584  track.set_py(total_momentum * dir.y());
585  track.set_pz(total_momentum * dir.z());
586 
587  // insert in map
588  m_track_map->insert(&track);
589 
590  if (Verbosity())
591  {
592  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
593  }
594  }
595 
596  // find collision point
597  /*
598  * intersection to either central membrane or endcaps
599  * if the position along beam and laser direction have the same sign, it will intercept the endcap
600  * otherwise will intercept the central membrane
601  */
602  const auto plane_strike = (pos.z() * dir.z() > 0) ? endcap_intersection(pos, dir) : central_membrane_intersection(pos, dir);
603 
604  // field cage intersection
605  const auto fc_strike = field_cage_intersection(pos, dir);
606 
607  // if none of the strikes is valid, there is no valid information found.
608  if (!(plane_strike || fc_strike)) return;
609 
610  // decide relevant end of laser
611  /* chose field cage intersection if valid, and if either plane intersection is invalid or happens on a larger z along the laser direction) */
612  const TVector3& strike = (fc_strike && (!plane_strike || fc_strike->z() / dir.z() < plane_strike->z() / dir.z())) ? *fc_strike : *plane_strike;
613 
614  //find length
615  TVector3 delta = (strike - pos);
616  double fullLength = delta.Mag();
617  int nHitSteps = fullLength / maxHitLength + 1;
618 
619  TVector3 start = pos;
620  TVector3 end = start;
621  TVector3 step = dir * (maxHitLength / (dir.Mag()));
622  double stepLength = 0;
623 
624  if (Verbosity())
625  {
626  std::cout << "PHG4TpcDirectLaser::AppendLaserTrack -"
627  << " fullLength: " << fullLength
628  << " nHitSteps: " << nHitSteps
629  << std::endl;
630  }
631 
632  for (int i = 0; i < nHitSteps; i++)
633  {
634  start = end; //new starting point is the previous ending point.
635  if (i + 1 == nHitSteps)
636  {
637  //last step is the remainder size
638  end = strike;
639  delta = start - end;
640  stepLength = delta.Mag();
641  }
642  else
643  {
644  //all other steps are uniform length
645  end = start + step;
646  stepLength = step.Mag();
647  }
648 
649  //from phg4tpcsteppingaction.cc
650  auto hit = new PHG4Hit_t;
651  hit->set_trkid(trackid);
652  hit->set_layer(99);
653 
654  //here we set the entrance values in cm
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);
659 
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);
664 
665  // momentum
666  hit->set_px(0, dir.X()); // GeV
667  hit->set_py(0, dir.Y());
668  hit->set_pz(0, dir.Z());
669 
670  hit->set_px(1, dir.X());
671  hit->set_py(1, dir.Y());
672  hit->set_pz(1, dir.Z());
673 
674  const double totalE = electrons_per_cm * stepLength / electrons_per_gev;
675 
676  hit->set_eion(totalE);
677  hit->set_edep(totalE);
678  m_g4hitcontainer->AddHit(detId, hit);
679  }
680 
681  return;
682 }