sPhenix code displayed by LXR master/​coresoftware/​simulation/​g4simulation/​g4tpc/​PHG4TpcDirectLaser.cc	Source navigation Diff markup Identifier search General search
	Version: master Revert   Change

File indexing completed on 2025-12-17 09:22:11

 #include "PHG4TpcDirectLaser.h"
 
 #include <phparameter/PHParameterInterface.h>  // for PHParameterInterface
 
 #include <g4main/PHG4HitContainer.h>
 #include <g4main/PHG4HitDefs.h>  // for get_volume_id
 #include <g4main/PHG4Hitv1.h>
 #include <g4main/PHG4Particlev3.h>
 #include <g4main/PHG4TruthInfoContainer.h>
 #include <g4main/PHG4VtxPointv1.h>
 
 #include <fun4all/Fun4AllReturnCodes.h>
 #include <fun4all/SubsysReco.h>  // for SubsysReco
 
 #include <phool/PHCompositeNode.h>
 #include <phool/PHIODataNode.h>
 #include <phool/PHNode.h>
 #include <phool/PHNodeIterator.h>
 #include <phool/PHObject.h>  // for PHObject
 #include <phool/getClass.h>
 #include <phool/phool.h>  // for PHWHERE
 
 #include <g4detectors/PHG4TpcGeom.h>
 #include <g4detectors/PHG4TpcGeomContainer.h>
 
 #include <trackbase_historic/SvtxTrackMap.h>
 #include <trackbase_historic/SvtxTrackMap_v2.h>
 #include <trackbase_historic/SvtxTrack_v2.h>
 
 #include <TFile.h>
 #include <TNtuple.h>
 #include <TVector3.h>  // for TVector3, operator*
 
 #include <gsl/gsl_const_mksa.h>  // for the speed of light
 
 #include <cassert>
 #include <iostream>  // for operator<<, basic_os...
 #include <optional>
 
 namespace
 {
   using PHG4Particle_t = PHG4Particlev3;
   using PHG4VtxPoint_t = PHG4VtxPointv1;
   using PHG4Hit_t = PHG4Hitv1;
 
   // utility
   template <class T>
   constexpr T square(const T& x)
   {
     return x * x;
   }
 
   // unique detector id for all direct lasers
   const int detId = PHG4HitDefs::get_volume_id("PHG4TpcDirectLaser");
 
   ///@name units
   //@{
   constexpr double cm = 1.0;
   //@}
 
   /// speed of light, in cm per ns
   constexpr double speed_of_light = GSL_CONST_MKSA_SPEED_OF_LIGHT * 1e-7;
 
   /// length of generated G4Hits along laser track
   constexpr double maxHitLength = 1. * cm;
 
   // inner and outer radii of field cages/TPC
   constexpr double begin_CM = 20. * cm;
   constexpr double end_CM = 78. * cm;
 
   //_____________________________________________________________
   std::optional<TVector3> central_membrane_intersection(const TVector3& start, const TVector3& direction, double halfwidth_CM)
   {
     const double end = start.z() > 0 ? halfwidth_CM : -halfwidth_CM;
     const double dist = end - start.z();
 
     // if line is vertical, it will never intercept the endcap
     if (direction.z() == 0)
     {
       return std::nullopt;
     }
 
     // check that distance and direction have the same sign
     if (dist * direction.z() < 0)
     {
       return std::nullopt;
     }
 
     const double direction_scale = dist / direction.z();
     return start + direction * direction_scale;
   }
 
   //_____________________________________________________________
   std::optional<TVector3> endcap_intersection(const TVector3& start, const TVector3& direction, double halflength_tpc)
   {
     const double end = start.z() > 0 ? halflength_tpc : -halflength_tpc;
     const double dist = end - start.z();
 
     // if line is vertical, it will never intercept the endcap
     if (direction.z() == 0)
     {
       return std::nullopt;
     }
 
     // check that distance and direction have the same sign
     if (dist * direction.z() < 0)
     {
       return std::nullopt;
     }
 
     const double direction_scale = dist / direction.z();
     return start + direction * direction_scale;
   }
 
   //_____________________________________________________________
   std::optional<TVector3> cylinder_line_intersection(const TVector3& s, const TVector3& v, double radius)
   {
     const double R2 = square(radius);
 
     // Generalized Parameters for collision with cylinder of radius R:
     // from quadratic formula solutions of when a vector intersects a circle:
     const double a = square(v.x()) + square(v.y());
     const double b = 2 * (v.x() * s.x() + v.y() * s.y());
     const double c = square(s.x()) + square(s.y()) - R2;
 
     const double rootterm = square(b) - 4 * a * c;
 
     /*
      * if a==0 then we are parallel and will have no solutions.
      * if the rootterm is negative, we will have no real roots,
      * we are outside the cylinder and pointing skew to the cylinder such that we never cross.
      */
     if (rootterm < 0 || a == 0)
     {
       return std::nullopt;
     }
 
     // Find the (up to) two points where we collide with the cylinder:
     const double sqrtterm = std::sqrt(rootterm);
     const double t1 = (-b + sqrtterm) / (2 * a);
     const double t2 = (-b - sqrtterm) / (2 * a);
 
     /*
      * if either of the t's are nonzero, we have a collision
      * the collision closest to the start (hence with the smallest t that is greater than zero) is the one that happens.
      */
     const double& min_t = (t2 < t1 && t2 > 0) ? t2 : t1;
     return s + v * min_t;
   }
 
   //_____________________________________________________________
   std::optional<TVector3> field_cage_intersection(const TVector3& start, const TVector3& direction)
   {
     auto ofc_strike = cylinder_line_intersection(start, direction, end_CM);
     auto ifc_strike = cylinder_line_intersection(start, direction, begin_CM);
 
     // if either of the two intersection is invalid, return the other
     if (!ifc_strike)
     {
       return ofc_strike;
     }
     if (!ofc_strike)
     {
       return ifc_strike;
     }
 
     // both intersection are valid, calculate signed distance to start z
     const auto ifc_dist = (ifc_strike->Z() - start.Z()) / direction.Z();
     const auto ofc_dist = (ofc_strike->Z() - start.Z()) / direction.Z();
 
     if (ifc_dist < 0)
     {
       return (ofc_dist > 0) ? ofc_strike : std::nullopt;
     }
     if (ofc_dist < 0)
     {
       return ifc_strike;
     }
 
     return (ifc_dist < ofc_dist) ? ifc_strike : ofc_strike;
   }
 
   /// TVector3 stream
   inline std::ostream& operator<<(std::ostream& out, const TVector3& vector)
   {
     out << "( " << vector.x() << ", " << vector.y() << ", " << vector.z() << ")";
     return out;
   }
 
 }  // namespace
 
 //_____________________________________________________________
 PHG4TpcDirectLaser::PHG4TpcDirectLaser(const std::string& name)
   : SubsysReco(name)
   , PHParameterInterface(name)
 {
   InitializeParameters();
 }
 
 //_____________________________________________________________
 int PHG4TpcDirectLaser::InitRun(PHCompositeNode* topNode)
 {
   // g4 truth info
   m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
   if (!m_g4truthinfo)
   {
     std::cout << "Fun4AllDstPileupMerger::load_nodes - creating node G4TruthInfo" << std::endl;
 
     PHNodeIterator iter(topNode);
     auto* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
     if (!dstNode)
     {
       std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
       return Fun4AllReturnCodes::ABORTRUN;
     }
 
     m_g4truthinfo = new PHG4TruthInfoContainer();
     dstNode->addNode(new PHIODataNode<PHObject>(m_g4truthinfo, "G4TruthInfo", "PHObject"));
   }
 
   // load and check G4Hit node
   hitnodename = "G4HIT_" + detector;
   auto* g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
   if (!g4hit)
   {
     std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
     return Fun4AllReturnCodes::ABORTRUN;
   }
 
   // find or create track map
   /* it is used to store laser parameters on a per event basis */
   m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
   if (!m_track_map)
   {
     // find DST node and check
     PHNodeIterator iter(topNode);
     auto* dstNode = static_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
     if (!dstNode)
     {
       std::cout << PHWHERE << "DST Node missing, aborting." << std::endl;
       return Fun4AllReturnCodes::ABORTRUN;
     }
 
     // find or create SVTX node
     iter = PHNodeIterator(dstNode);
     auto* node = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "SVTX"));
     if (!node)
     {
       dstNode->addNode(node = new PHCompositeNode("SVTX"));
     }
 
     // add track node
     m_track_map = new SvtxTrackMap_v2;
     node->addNode(new PHIODataNode<PHObject>(m_track_map, m_track_map_name, "PHObject"));
   }
 
   // setup parameters
   UpdateParametersWithMacro();
   electrons_per_cm = get_int_param("electrons_per_cm");
   electrons_per_gev = get_double_param("electrons_per_gev");
 
     PHG4TpcGeomContainer *GeomContainer = findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
    
   PHG4TpcGeom *layergeom = GeomContainer->GetLayerCellGeom(20);  // z geometry is the same for all layers
   double maxdriftlength = layergeom->get_max_driftlength();
   halfwidth_CM = layergeom->get_CM_halfwidth();
   halflength_tpc = maxdriftlength + halfwidth_CM;
   
   // setup lasers
   SetupLasers();
 
   // print configuration
   if (m_steppingpattern == true)
   {
     std::cout << "PHG4TpcDirectLaser::InitRun - m_steppingpattern: " << m_steppingpattern << std::endl;
     std::cout << "PHG4TpcDirectLaser::InitRun - nTotalSteps: " << nTotalSteps << std::endl;
   }
   else
   {
     std::cout << "PHG4TpcDirectLaser::InitRun - m_autoAdvanceDirectLaser: " << m_autoAdvanceDirectLaser << std::endl;
     std::cout << "PHG4TpcDirectLaser::InitRun - phi steps: " << nPhiSteps << " min: " << minPhi << " max: " << maxPhi << std::endl;
     std::cout << "PHG4TpcDirectLaser::InitRun - theta steps: " << nThetaSteps << " min: " << minTheta << " max: " << maxTheta << std::endl;
     std::cout << "PHG4TpcDirectLaser::InitRun - nTotalSteps: " << nTotalSteps << std::endl;
   }
   std::cout << "PHG4TpcDirectLaser::InitRun - electrons_per_cm: " << electrons_per_cm << std::endl;
   std::cout << "PHG4TpcDirectLaser::InitRun - electrons_per_gev " << electrons_per_gev << std::endl;
 
   // TFile * infile1 = TFile::Open("theta_phi_laser.root");
 
   std::string LASER_ANGLES_ROOTFILE = std::string(getenv("CALIBRATIONROOT")) + "/TPC/DirectLaser/theta_phi_laser.root";
   TFile* infile1 = TFile::Open(LASER_ANGLES_ROOTFILE.c_str());
 
   pattern = (TNtuple*) infile1->Get("angles");
   pattern->SetBranchAddress("#theta", &theta_p);
   pattern->SetBranchAddress("#phi", &phi_p);
 
   return Fun4AllReturnCodes::EVENT_OK;
 }
 
 //_____________________________________________________________
 int PHG4TpcDirectLaser::process_event(PHCompositeNode* topNode)
 {
   // g4 input event
   m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
   assert(m_g4truthinfo);
 
   // load g4hit container
   m_g4hitcontainer = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
   assert(m_g4hitcontainer);
 
   // load track map
   m_track_map = findNode::getClass<SvtxTrackMap>(topNode, m_track_map_name);
   assert(m_track_map);
 
   if (m_autoAdvanceDirectLaser || m_steppingpattern)
   {
     AimToNextPatternStep();
   }
   //_________________________________________________
   else
   {
     // use arbitrary direction
     AimToThetaPhi(arbitrary_theta, arbitrary_phi);
   }
 
   return Fun4AllReturnCodes::EVENT_OK;
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::SetDefaultParameters()
 {
   // same gas parameters as in PHG4TpcElectronDrift::SetDefaultParameters
 
   // Data on gasses @20 C and 760 Torr from the following source:
   // http://www.slac.stanford.edu/pubs/icfa/summer98/paper3/paper3.pdf
   // diffusion and drift velocity for 400kV for NeCF4 50/50 from calculations:
   // http://skipper.physics.sunysb.edu/~prakhar/tpc/HTML_Gases/split.html
   static constexpr double Ne_dEdx = 1.56;    // keV/cm
   static constexpr double CF4_dEdx = 7.00;   // keV/cm
   static constexpr double Ne_NTotal = 43;    // Number/cm
   static constexpr double CF4_NTotal = 100;  // Number/cm
   static constexpr double Tpc_NTot = 0.5 * Ne_NTotal + 0.5 * CF4_NTotal;
   static constexpr double Tpc_dEdx = 0.5 * Ne_dEdx + 0.5 * CF4_dEdx;
   static constexpr double Tpc_ElectronsPerKeV = Tpc_NTot / Tpc_dEdx;
 
   // number of electrons per deposited GeV in TPC gas
   set_default_double_param("electrons_per_gev", Tpc_ElectronsPerKeV * 1e6);
 
   //  set_default_double_param("tpc_half_length", 102.325);
   //  set_default_double_param("CM_halfwidth", 0.28);
 
   // number of electrons deposited by laser per cm
   set_default_int_param("electrons_per_cm", 72);
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::SetPhiStepping(int n, double min, double max)
 {
   if (n < 0 || max < min)
   {
     std::cout << PHWHERE << " - invalid" << std::endl;
     return;
   }
   nPhiSteps = n;
   minPhi = min;
   maxPhi = max;
   nTotalSteps = nThetaSteps * nPhiSteps;
   return;
 }
 //_____________________________________________________________
 void PHG4TpcDirectLaser::SetThetaStepping(int n, double min, double max)
 {
   if (n < 0 || max < min)
   {
     std::cout << PHWHERE << " - invalid" << std::endl;
     return;
   }
   nThetaSteps = n;
   minTheta = min;
   maxTheta = max;
   nTotalSteps = nThetaSteps * nPhiSteps;
 
   return;
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::SetFileStepping(int n)
 {
   if (n < 0 || n > 13802)  // 13802 = hard coded number of tuple entries
   {
     std::cout << PHWHERE << " - invalid" << std::endl;
     return;
   }
   nTotalSteps = n;
 
   return;
 }
 
 //_____________________________________________________________
 
 void PHG4TpcDirectLaser::SetupLasers()
 {
   // clear previous lasers
   m_lasers.clear();
 
   // position of first laser at positive z
   const TVector3 position_base(60 * cm, 0., halflength_tpc);
 
   // add lasers
   for (int i = 0; i < 8; ++i)
   {
     Laser laser;
 
     // set laser direction
     /*
      * first four lasers are on positive z readout plane, and shoot towards negative z
      * next four lasers are on negative z readout plane and shoot towards positive z
      */
     laser.m_position = position_base;
     if (i < 4)
     {
       laser.m_position.SetZ(position_base.z());
       laser.m_direction = -1;
       laser.m_phi = M_PI / 2 * i - (15 * M_PI / 180);  // additional offset of 15 deg.
     }
     else
     {
       laser.m_position.SetZ(-position_base.z());
       laser.m_direction = 1;
       laser.m_phi = M_PI / 2 * i + (15 * M_PI / 180);  // additional offset of 15 deg.
     }
 
     // rotate around z
     laser.m_position.RotateZ(laser.m_phi);
 
     // append
     m_lasers.push_back(laser);  // All lasers
     //  if(i==0) m_lasers.push_back(laser);//Only laser 1
     //  if(i==3) m_lasers.push_back(laser);// Laser 4
     // if(i<4) m_lasers.push_back(laser);//Lasers 1, 2, 3, 4
   }
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::AimToNextPatternStep()
 {
   if (nTotalSteps >= 1)
   {
     if (m_steppingpattern)
     {
       AimToPatternStep_File(currentPatternStep);
       ++currentPatternStep;
     }
     else
     {
       AimToPatternStep(currentPatternStep);
       ++currentPatternStep;
     }
   }
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::AimToThetaPhi(double theta, double phi)
 {
   if (Verbosity())
   {
     std::cout << "PHG4TpcDirectLaser::AimToThetaPhi - theta: " << theta << " phi: " << phi << std::endl;
   }
 
   // all lasers
   for (const auto& laser : m_lasers)
   {
     AppendLaserTrack(theta, phi, laser);
   }
 }
 
 //_____________________________________________________________
 void PHG4TpcDirectLaser::AimToPatternStep(int n)
 {
   // trim against overflows
   n = n % nTotalSteps;
 
   if (Verbosity())
   {
     std::cout << "PHG4TpcDirectLaser::AimToPatternStep - step: " << n << "/" << nTotalSteps << std::endl;
   }
 
   // store as current pattern
   currentPatternStep = n;
 
   // calculate theta
   const int thetaStep = n / nPhiSteps;
   const double theta = minTheta + thetaStep * (maxTheta - minTheta) / nThetaSteps;
 
   // calculate phi
   const int phiStep = n % nPhiSteps;
   const double phi = minPhi + phiStep * (maxPhi - minPhi) / nPhiSteps;
 
   // generate laser tracks
   AimToThetaPhi(theta, phi);
 
   return;
 }
 
 //_____________________________________________________________
 
 void PHG4TpcDirectLaser::AimToPatternStep_File(int n)
 {
   // trim against overflows
   n = n % nTotalSteps;
 
   if (Verbosity())
   {
     std::cout << "PHG4TpcDirectLaser::AimToPatternStep_File - step: " << n << "/" << nTotalSteps << std::endl;
   }
 
   // store as current pattern
   currentPatternStep = n;
 
   pattern->GetEntry(n);
 
   // calculate theta
   std::cout << "From file, current entry = " << n << " Theta: " << theta_p << " Phi: " << phi_p << std::endl;
 
   const double theta = theta_p * M_PI / 180.;
 
   // calculate phi
   const double phi = phi_p * M_PI / 180.;
 
   // generate laser tracks
   AimToThetaPhi(theta, phi);
 
   return;
 }
 
 //_____________________________________________________________
 
 void PHG4TpcDirectLaser::AppendLaserTrack(double theta, double phi, const PHG4TpcDirectLaser::Laser& laser)
 {
   if (!m_g4hitcontainer)
   {
     std::cout << PHWHERE << "invalid g4hit container. aborting" << std::endl;
     return;
   }
 
   // store laser position
   const auto& pos = laser.m_position;
 
   // define track direction
   const auto& direction = laser.m_direction;
   TVector3 dir(0, 0, direction);
 
   // adjust direction
   dir.RotateY(theta * direction);
 
   if (laser.m_direction == -1)
   {
     dir.RotateZ(phi);  // if +z facing -z
   }
   else
   {
     dir.RotateZ(-phi);  // if -z facting +z
   }
 
   // also rotate by laser azimuth
   dir.RotateZ(laser.m_phi);
 
   // print
   if (Verbosity())
   {
     std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
   }
 
   // dummy momentum
   static constexpr double total_momentum = 1;
 
   // mc track id
   int trackid = -1;
 
   // create truth vertex and particle
   if (m_g4truthinfo)
   {
     // add vertex
     const auto vtxid = m_g4truthinfo->maxvtxindex() + 1;
     auto* const vertex = new PHG4VtxPoint_t(pos.x(), pos.y(), pos.z(), 0, vtxid);
     m_g4truthinfo->AddVertex(vtxid, vertex);
 
     // increment track id
     trackid = m_g4truthinfo->maxtrkindex() + 1;
 
     // create new g4particle
     auto* particle = new PHG4Particle_t();
     particle->set_track_id(trackid);
     particle->set_vtx_id(vtxid);
     particle->set_parent_id(0);
     particle->set_primary_id(trackid);
     particle->set_px(total_momentum * dir.x());
     particle->set_py(total_momentum * dir.y());
     particle->set_pz(total_momentum * dir.z());
 
     m_g4truthinfo->AddParticle(trackid, particle);
   }
 
   // store in SvtxTrack map
   if (m_track_map)
   {
     SvtxTrack_v2 track;
     track.set_x(pos.x());
     track.set_y(pos.y());
     track.set_z(pos.z());
 
     // total momentum is irrelevant. What matters is the direction
     track.set_px(total_momentum * dir.x());
     track.set_py(total_momentum * dir.y());
     track.set_pz(total_momentum * dir.z());
 
     // insert in map
     m_track_map->insert(&track);
 
     if (Verbosity())
     {
       std::cout << "PHG4TpcDirectLaser::AppendLaserTrack - position: " << pos << " direction: " << dir << std::endl;
     }
   }
 
   // find collision point
   /*
    * intersection to either central membrane or endcaps
    * if the position along beam and laser direction have the same sign, it will intercept the endcap
    * otherwise will intercept the central membrane
    */
   const auto plane_strike = (pos.z() * dir.z() > 0) ? endcap_intersection(pos, dir,halflength_tpc) : central_membrane_intersection(pos, dir, halfwidth_CM);
 
   // field cage intersection
   const auto fc_strike = field_cage_intersection(pos, dir);
 
   // if none of the strikes is valid, there is no valid information found.
   if (!(plane_strike || fc_strike))
   {
     return;
   }
 
   // decide relevant end of laser
   /* chose field cage intersection if valid, and if either plane intersection is invalid or happens on a larger z along the laser direction) */
   const TVector3& strike = (fc_strike && (!plane_strike || fc_strike->z() / dir.z() < plane_strike->z() / dir.z())) ? *fc_strike : *plane_strike;
 
   // find length
   TVector3 delta = (strike - pos);
   double fullLength = delta.Mag();
   int nHitSteps = fullLength / maxHitLength + 1;
 
   TVector3 start = pos;
   TVector3 end = start;
   TVector3 step = dir * (maxHitLength / (dir.Mag()));
   double stepLength = 0;
 
   if (Verbosity())
   {
     std::cout << "PHG4TpcDirectLaser::AppendLaserTrack -"
               << " fullLength: " << fullLength
               << " nHitSteps: " << nHitSteps
               << std::endl;
   }
 
   for (int i = 0; i < nHitSteps; i++)
   {
     start = end;  // new starting point is the previous ending point.
     if (i + 1 == nHitSteps)
     {
       // last step is the remainder size
       end = strike;
       delta = start - end;
       stepLength = delta.Mag();
     }
     else
     {
       // all other steps are uniform length
       end = start + step;
       stepLength = step.Mag();
     }
 
     // from phg4tpcsteppingaction.cc
     auto* hit = new PHG4Hit_t;
     hit->set_trkid(trackid);
     hit->set_layer(99);
 
     // here we set the entrance values in cm
     hit->set_x(0, start.X() / cm);
     hit->set_y(0, start.Y() / cm);
     hit->set_z(0, start.Z() / cm);
     hit->set_t(0, (start - pos).Mag() / speed_of_light);
 
     hit->set_x(1, end.X() / cm);
     hit->set_y(1, end.Y() / cm);
     hit->set_z(1, end.Z() / cm);
     hit->set_t(1, (end - pos).Mag() / speed_of_light);
 
     // momentum
     hit->set_px(0, dir.X());  // GeV
     hit->set_py(0, dir.Y());
     hit->set_pz(0, dir.Z());
 
     hit->set_px(1, dir.X());
     hit->set_py(1, dir.Y());
     hit->set_pz(1, dir.Z());
 
     const double totalE = electrons_per_cm * stepLength / electrons_per_gev;
 
     hit->set_eion(totalE);
     hit->set_edep(totalE);
     m_g4hitcontainer->AddHit(detId, hit);
   }
 
   return;
 }

This page was automatically generated by the 2.3.7 LXR engine. The LXR team