Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:21 

0001 
0002 #include "TrackSeedTrackMapConverter.h"
0003 
0004 #include <trackbase/ActsGeometry.h>
0005 #include <trackbase/TrackFitUtils.h>
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrClusterContainer.h>
0008 
0009 #include <trackbase_historic/SvtxTrackMap_v1.h>
0010 #include <trackbase_historic/SvtxTrack_v4.h>
0011 #include <trackbase_historic/SvtxTrackState_v2.h>
0012 #include <trackbase_historic/TrackSeed.h>
0013 #include <trackbase_historic/TrackSeedContainer.h>
0014 #include <trackbase_historic/TrackSeedHelper.h>
0015 
0016 #include <fun4all/Fun4AllReturnCodes.h>
0017 
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>  // for PHNode
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>  // for PHObject
0023 #include <phool/getClass.h>
0024 #include <phool/phool.h>
0025 
0026 #include <cmath>
0027 namespace
0028 {
0029   template <class T>
0030   inline T square(const T& x)
0031   {
0032     return x * x;
0033   }
0034 }  // namespace
0035 //____________________________________________________________________________..
0036 TrackSeedTrackMapConverter::TrackSeedTrackMapConverter(const std::string& name)
0037   : SubsysReco(name)
0038 {
0039 }
0040 
0041 //____________________________________________________________________________..
0042 int TrackSeedTrackMapConverter::InitRun(PHCompositeNode* topNode)
0043 {
0044   int ret = getNodes(topNode);
0045   std::istringstream stringline(m_fieldMap);
0046   stringline >> fieldstrength;
0047   if (!stringline.fail())  // it is a float
0048   {
0049     m_ConstField = true;
0050     if(std::stof(m_fieldMap) < 0.02)
0051     {
0052       m_zeroField = true;
0053     }
0054   }
0055   return ret;
0056 }
0057 
0058 //____________________________________________________________________________..
0059 int TrackSeedTrackMapConverter::process_event(PHCompositeNode* /*unused*/)
0060 {
0061   if (Verbosity() > 1)
0062   {
0063     std::cout << "silicon seed map size " << m_siContainer->size() << std::endl;
0064     for (auto iter = m_siContainer->begin(); iter != m_siContainer->end();
0065          ++iter)
0066     {
0067       auto seed = *iter;
0068       if (!seed)
0069       {
0070         std::cout << "no seed at index " << m_siContainer->index(iter)
0071                   << std::endl;
0072         continue;
0073       }
0074       seed->identify();
0075     }
0076     if (m_tpcContainer)
0077     {
0078       std::cout << "tpc seed map size " << m_tpcContainer->size() << std::endl;
0079       for (auto iter = m_tpcContainer->begin(); iter != m_tpcContainer->end();
0080            ++iter)
0081       {
0082         auto seed = *iter;
0083         if (!seed)
0084         {
0085           std::cout << "no tpc seed at entry " << m_tpcContainer->index(iter)
0086                     << std::endl;
0087           continue;
0088         }
0089         seed->identify();
0090       }
0091     }
0092   }
0093 
0094   /// Start with a fresh track map in case running iterative tracking
0095   m_trackMap->Reset();
0096 
0097   unsigned int trackid = 0;
0098   for (const auto& trackSeed : *m_seedContainer)
0099   {
0100     /// If the seed was removed, skip it
0101     if (!trackSeed)
0102     {
0103       continue;
0104     }
0105 
0106     if (m_trackSeedName.find("SvtxTrackSeed") != std::string::npos)
0107     {
0108       /// Catches entries in the vector removed by ghost finder
0109       unsigned int tpcindex = trackSeed->get_tpc_seed_index();
0110       TrackSeed* seed = m_tpcContainer->get(tpcindex);
0111       if (!seed)
0112       {
0113         continue;
0114       }
0115     }
0116 
0117     auto svtxtrack = std::make_unique<SvtxTrack_v4>();
0118 
0119     if (Verbosity() > 0)
0120     {
0121       std::cout << "iterating track seed " << trackid << std::endl;
0122     }
0123 
0124     svtxtrack->set_id(trackid);
0125     trackid++;
0126     /// If we've run the track matching
0127     if (m_trackSeedName.find("SvtxTrackSeed") != std::string::npos)
0128     {
0129       if (Verbosity() > 0)
0130       {
0131         std::cout << "tpc seed id " << trackSeed->get_tpc_seed_index() << std::endl;
0132         std::cout << "si seed id " << trackSeed->get_silicon_seed_index() << std::endl;
0133         std::cout << "crossing_estimate " << trackSeed->get_crossing_estimate() << std::endl;
0134       }
0135 
0136       unsigned int seedindex = trackSeed->get_tpc_seed_index();
0137       TrackSeed* tpcseed = m_tpcContainer->get(seedindex);
0138       if (!m_cosmics)
0139       {
0140         if (trackSeed->get_silicon_seed_index() == std::numeric_limits<unsigned int>::max())
0141         {
0142           /// Didn't find a match, so just use the tpc seed
0143           const auto position = TrackSeedHelper::get_xyz(tpcseed);
0144           svtxtrack->set_x(position.x());
0145           svtxtrack->set_y(position.y());
0146           svtxtrack->set_z(position.z());
0147           svtxtrack->set_crossing(0);
0148         }
0149         else
0150         {
0151           TrackSeed* siseed = m_siContainer->get(trackSeed->get_silicon_seed_index());
0152           const auto position = TrackSeedHelper::get_xyz(siseed);
0153           svtxtrack->set_x(position.x());
0154           svtxtrack->set_y(position.y());
0155           svtxtrack->set_z(position.z());
0156           svtxtrack->set_crossing(siseed->get_crossing());
0157 
0158           svtxtrack->set_silicon_seed(siseed);
0159         }
0160         svtxtrack->set_charge(tpcseed->get_qOverR() > 0 ? 1 : -1);
0161         if (!m_ConstField)
0162         {
0163           svtxtrack->set_px(tpcseed->get_pt() * std::cos(tpcseed->get_phi()));
0164           svtxtrack->set_py(tpcseed->get_pt() * std::sin(tpcseed->get_phi()));
0165           svtxtrack->set_pz(tpcseed->get_pz());
0166         }
0167         else
0168         {
0169           float pt = fabs(1. / tpcseed->get_qOverR()) * (0.3 / 100) * fieldstrength;
0170           float phi = tpcseed->get_phi();
0171           svtxtrack->set_px(pt * std::cos(phi));
0172           svtxtrack->set_py(pt * std::sin(phi));
0173           svtxtrack->set_pz(pt * std::cosh(tpcseed->get_eta()) * std::cos(tpcseed->get_theta()));
0174         }
0175       }
0176       else
0177       {
0178         unsigned int silseedindex = trackSeed->get_silicon_seed_index();
0179         TrackSeed* silseed = m_siContainer->get(silseedindex);
0180 
0181         // get positions from cluster keys
0182         // TODO: should implement distortions
0183         TrackSeedHelper::position_map_t positions;
0184         for( auto key_iter = tpcseed->begin_cluster_keys(); key_iter != tpcseed->end_cluster_keys(); ++key_iter )
0185         {
0186           const auto& key(*key_iter);
0187           positions.emplace(key, m_tGeometry->getGlobalPosition( key, m_clusters->findCluster(key)));
0188         }
0189 
0190         // perform circle fit
0191         TrackSeedHelper::circleFitByTaubin(tpcseed, positions, 0, 58);
0192 
0193         float tpcR = fabs(1. / tpcseed->get_qOverR());
0194         float tpcx = tpcseed->get_X0();
0195         float tpcy = tpcseed->get_Y0();
0196         float vertexradius = 80;
0197         const auto intersect = TrackFitUtils::circle_circle_intersection(vertexradius, tpcR, tpcx, tpcy);
0198         float intx, inty;
0199 
0200         if (std::get<1>(intersect) < std::get<3>(intersect))
0201         {
0202           intx = std::get<0>(intersect);
0203           inty = std::get<1>(intersect);
0204         }
0205         else
0206         {
0207           intx = std::get<2>(intersect);
0208           inty = std::get<3>(intersect);
0209         }
0210         float slope = tpcseed->get_slope();
0211 
0212         float intz = vertexradius * slope + tpcseed->get_Z0();
0213 
0214         Acts::Vector3 inter(intx, inty, intz);
0215 
0216         std::vector<float> tpcparams{tpcR, tpcx, tpcy, tpcseed->get_slope(),
0217                                      tpcseed->get_Z0()};
0218         auto tangent = TrackFitUtils::get_helix_tangent(tpcparams,
0219                                                         inter);
0220 
0221         auto tan = tangent.second;
0222         auto pca = tangent.first;
0223 
0224         svtxtrack->set_x(pca.x());
0225         svtxtrack->set_y(pca.y());
0226 
0227         float p;
0228         if (m_ConstField)
0229         {
0230           p = std::cosh(tpcseed->get_eta()) * fabs(1. / tpcseed->get_qOverR()) * (0.3 / 100) * fieldstrength;
0231         }
0232         else
0233         {
0234           p = tpcseed->get_p();
0235         }
0236 
0237         tan *= p;
0238         auto [charge, cosmicslope] = getCosmicCharge(tpcseed, vertexradius);
0239 
0240         svtxtrack->set_px(charge < 0 ? tan.x() : tan.x() * -1);
0241         svtxtrack->set_py(charge < 0 ? tan.y() : tan.y() * -1);
0242         svtxtrack->set_pz(cosmicslope > 0 ? fabs(tan.z()) : -1 * fabs(tan.z()));
0243         svtxtrack->set_z(svtxtrack->get_pz() > 0 ? (slope < 0 ? intz : vertexradius * slope * -1 + tpcseed->get_Z0()) : (slope > 0 ? intz : vertexradius * slope * -1 + tpcseed->get_Z0()));
0244         svtxtrack->set_charge(charge);
0245         addKeys(svtxtrack, tpcseed);
0246         if (silseed)
0247         {
0248           addKeys(svtxtrack, silseed);
0249         }
0250 
0251         svtxtrack->set_tpc_seed(tpcseed);
0252         svtxtrack->set_silicon_seed(silseed);
0253         if (Verbosity() > 5)
0254         {
0255           svtxtrack->identify();
0256         }
0257       }
0258       svtxtrack->set_tpc_seed(tpcseed);
0259     }
0260 
0261     else
0262     {
0263       /// Otherwise we are using an individual subdetectors container
0264       svtxtrack->set_id(m_seedContainer->find(trackSeed));
0265 
0266       const auto position = TrackSeedHelper::get_xyz(trackSeed);
0267       svtxtrack->set_x(position.x());
0268       svtxtrack->set_y(position.y());
0269       if(m_zeroField)
0270       {
0271         // replace x and y with a line fit instead of helix fit
0272         lineFit(svtxtrack.get(), trackSeed);
0273       }
0274       svtxtrack->set_z(position.z());
0275       svtxtrack->set_charge(trackSeed->get_qOverR() > 0 ? 1 : -1);
0276       if (m_ConstField)
0277       {
0278         float pt = fabs(1. / trackSeed->get_qOverR()) * (0.3 / 100) * fieldstrength;
0279 
0280         float phi = trackSeed->get_phi();
0281         svtxtrack->set_px(pt * std::cos(phi));
0282         svtxtrack->set_py(pt * std::sin(phi));
0283         svtxtrack->set_pz(pt * std::cosh(trackSeed->get_eta()) * std::cos(trackSeed->get_theta()));
0284       }
0285       else
0286       {
0287 
0288         svtxtrack->set_px(trackSeed->get_pt() * std::cos(trackSeed->get_phi()));
0289         svtxtrack->set_py(trackSeed->get_pt() * std::sin(trackSeed->get_phi()));
0290         svtxtrack->set_pz(trackSeed->get_pz());
0291       }
0292 
0293       // calculate chisq and ndf
0294       float R = 1. / std::fabs(trackSeed->get_qOverR());
0295       float X0 = trackSeed->get_X0();
0296       float Y0 = trackSeed->get_Y0();
0297       float Z0 = trackSeed->get_Z0();
0298       float slope = trackSeed->get_slope();
0299 
0300       std::vector<float> fitpars = {R, X0, Y0, slope, Z0};
0301 
0302       svtxtrack->set_crossing(trackSeed->get_crossing());
0303 
0304       SvtxTrackState_v2 initial_state(0.);
0305       initial_state.set_x(svtxtrack->get_x());
0306       initial_state.set_y(svtxtrack->get_y());
0307       initial_state.set_z(svtxtrack->get_z());
0308       initial_state.set_px(svtxtrack->get_px());
0309       initial_state.set_py(svtxtrack->get_py());
0310       initial_state.set_pz(svtxtrack->get_pz());
0311       svtxtrack->insert_state(&initial_state);
0312 
0313       // TPC-only tracks don't have crossing info, so we set crossing to 0
0314       // Silicon-only tracks may have crossing info, so a crossing of SHRT_MAX indicates a genuine ambiguity and we leave it as is
0315       if (m_trackSeedName.find("TpcTrackSeed") != std::string::npos)
0316       {
0317         if(svtxtrack->get_crossing() == SHRT_MAX)
0318         {
0319           svtxtrack->set_crossing(0);
0320         }
0321       }
0322 
0323       std::vector<double> xy_error2;
0324       std::vector<double> rz_error2;
0325       std::vector<double> xy_residuals;
0326       std::vector<double> rz_residuals;
0327       std::vector<double> x_circle;
0328       std::vector<double> y_circle;
0329       std::vector<double> z_line;
0330 
0331       // TrackFitUtils::get_helix_tangent sometimes returns a tangent vector in the opposite direction that we would want
0332       // To guard against this, we keep track of the previous momentum vector and un-flip the current momentum vector if it flips
0333       float last_px = initial_state.get_px();
0334       float last_py = initial_state.get_py();
0335       float last_pz = initial_state.get_pz();
0336 
0337       for (auto c_iter = trackSeed->begin_cluster_keys();
0338            c_iter != trackSeed->end_cluster_keys();
0339            ++c_iter)
0340       {
0341         TrkrDefs::cluskey key = *c_iter;
0342         TrkrCluster* c = m_clusters->findCluster(key);
0343         Acts::Vector3 pos = m_tGeometry->getGlobalPosition(key, c);
0344 
0345         auto surf = m_tGeometry->maps().getSurface(key,c);
0346         Acts::Vector3 surface_center = surf->center(m_tGeometry->geometry().getGeoContext()) * 0.1; // to cm
0347         Acts::Vector3 intersection = TrackFitUtils::get_helix_surface_intersection(surf,fitpars,pos,m_tGeometry);
0348         float pathlength = TrackFitUtils::get_helix_pathlength(fitpars,position,intersection);
0349         std::pair<Acts::Vector3,Acts::Vector3> tangent = TrackFitUtils::get_helix_tangent(fitpars,surface_center);
0350 
0351         SvtxTrackState_v2 state(pathlength);
0352         state.set_x(intersection.x());
0353         state.set_y(intersection.y());
0354         state.set_z(intersection.z());
0355         float p = trackSeed->get_p();
0356         float new_px = p*tangent.second.x();
0357         float new_py = p*tangent.second.y();
0358         float new_pz = p*tangent.second.z();
0359         if((last_px*new_px+last_py*new_py+last_pz*new_pz)<0.)
0360         {
0361           new_px = -new_px;
0362           new_py = -new_py;
0363           new_pz = -new_pz;
0364         }
0365         state.set_px(new_px);
0366         state.set_py(new_py);
0367         state.set_pz(new_pz);
0368         state.set_cluskey(key);
0369         svtxtrack->insert_state(&state);
0370 
0371         last_px = new_px;
0372         last_py = new_py;
0373         last_pz = new_pz;
0374 
0375         double x = pos(0);
0376         double y = pos(1);
0377         double z = pos(2);
0378         double r = sqrt(x*x+y*y);
0379         double dx = x - intersection.x();
0380         double dy = y - intersection.y();
0381         double dz = z - intersection.z();
0382         double dr = r - sqrt(intersection.x()*intersection.x()+intersection.y()*intersection.y());
0383         xy_residuals.push_back(sqrt(dx*dx+dy*dy));
0384         rz_residuals.push_back(sqrt(dr*dr+dz*dz));
0385 
0386         // ignoring covariance for simplicity
0387         xy_error2.push_back(c->getRPhiError()*c->getRPhiError());
0388         rz_error2.push_back(c->getZError()*c->getZError());
0389         double phi = atan2(dy, dx);
0390         x_circle.push_back(R * cos(phi) + X0);
0391         y_circle.push_back(R * sin(phi) + Y0);
0392         z_line.push_back(R * slope + Z0);
0393       }
0394       double chi2 = 0.;
0395       for (unsigned int i = 0; i < xy_residuals.size(); i++)
0396       {
0397         // method lifted from GPUTPCTrackParam::Filter
0398         chi2 += xy_residuals[i] * xy_residuals[i] / xy_error2[i] + rz_residuals[i] * rz_residuals[i] / rz_error2[i];
0399       }
0400       svtxtrack->set_chisq(chi2);
0401       // GPUTPCTrackParam initially sets NDF to -3 on first cluster and increments by 2 with every application of filter
0402       svtxtrack->set_ndf(2 * xy_residuals.size() - 5);
0403 
0404       addKeys(svtxtrack, trackSeed);
0405       if (m_trackSeedName.find("SiliconTrackSeed") != std::string::npos)
0406       {
0407         svtxtrack->set_silicon_seed(trackSeed);
0408         svtxtrack->set_tpc_seed(nullptr);
0409       }
0410       else if (m_trackSeedName.find("TpcTrackSeed") != std::string::npos)
0411       {
0412         svtxtrack->set_tpc_seed(trackSeed);
0413         svtxtrack->set_silicon_seed(nullptr);
0414       }
0415     }
0416 
0417     if (Verbosity() > 0)
0418     {
0419       std::cout << "Inserting svtxtrack into map " << std::endl;
0420       svtxtrack->identify();
0421     }
0422 
0423     m_trackMap->insert(svtxtrack.get());
0424   }
0425 
0426   return Fun4AllReturnCodes::EVENT_OK;
0427 }
0428 void TrackSeedTrackMapConverter::lineFit(SvtxTrack_v4* track, TrackSeed* seed)
0429 {
0430   auto firstclusterckey = *(seed->begin_cluster_keys());
0431   auto cluster = m_clusters->findCluster(firstclusterckey);
0432   auto glob = m_tGeometry->getGlobalPosition(firstclusterckey, cluster);
0433   float phi = seed->get_phi();
0434   float theta = seed->get_theta();
0435   Acts::Vector3 mom(std::cos(phi) * std::sin(theta),
0436                     std::sin(phi) * std::sin(theta), std::cos(theta));
0437 
0438   auto pca = TrackFitUtils::getPCALinePoint(Acts::Vector3::Zero(), mom, glob);
0439   track->set_x(pca.x());
0440   track->set_y(pca.y());
0441 }
0442 //____________________________________________________________________________..
0443 int TrackSeedTrackMapConverter::End(PHCompositeNode* /*unused*/)
0444 {
0445   return Fun4AllReturnCodes::EVENT_OK;
0446 }
0447 void TrackSeedTrackMapConverter::addKeys(TrackSeed* seedToAddTo,
0448                                          TrackSeed* seedToAdd)
0449 {
0450   for (auto iter = seedToAdd->begin_cluster_keys();
0451        iter != seedToAdd->end_cluster_keys();
0452        ++iter)
0453   {
0454     seedToAddTo->insert_cluster_key(*iter);
0455   }
0456 }
0457 void TrackSeedTrackMapConverter::addKeys(std::unique_ptr<SvtxTrack_v4>& track, TrackSeed* seed)
0458 {
0459   for (TrackSeed::ConstClusterKeyIter iter = seed->begin_cluster_keys();
0460        iter != seed->end_cluster_keys();
0461        ++iter)
0462   {
0463     track->insert_cluster_key(*iter);
0464   }
0465 }
0466 
0467 int TrackSeedTrackMapConverter::getNodes(PHCompositeNode* topNode)
0468 {
0469   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0470   if (!m_trackMap)
0471   {
0472     // create it
0473     PHNodeIterator iter(topNode);
0474 
0475     PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0476         "PHCompositeNode", "DST"));
0477     if (!dstNode)
0478     {
0479       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0480       return Fun4AllReturnCodes::ABORTEVENT;
0481     }
0482     PHNodeIterator iter_dst(dstNode);
0483 
0484     // Create the SVTX node
0485     PHCompositeNode* tb_node =
0486         dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode",
0487                                                           "SVTX"));
0488     if (!tb_node)
0489     {
0490       tb_node = new PHCompositeNode("SVTX");
0491       dstNode->addNode(tb_node);
0492       if (Verbosity() > 0)
0493       {
0494         std::cout << PHWHERE << "SVTX node added" << std::endl;
0495       }
0496     }
0497     m_trackMap = new SvtxTrackMap_v1;
0498     PHIODataNode<PHObject>* tracks_node =
0499         new PHIODataNode<PHObject>(m_trackMap, "SvtxTrackMap", "PHObject");
0500     tb_node->addNode(tracks_node);
0501   }
0502 
0503   if(m_trackMapName != "SvtxTrackMap")
0504   {
0505   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_trackMapName);
0506   if (!m_trackMap)
0507   {
0508     PHNodeIterator iter(topNode);
0509 
0510     PHCompositeNode* dstNode = static_cast<PHCompositeNode*>(iter.findFirst(
0511         "PHCompositeNode", "DST"));
0512     if (!dstNode)
0513     {
0514       std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0515       return Fun4AllReturnCodes::ABORTEVENT;
0516     }
0517     PHNodeIterator iter_dst(dstNode);
0518 
0519     // Create the SVTX node
0520     PHCompositeNode* tb_node =
0521         dynamic_cast<PHCompositeNode*>(iter_dst.findFirst("PHCompositeNode",
0522                                                           "SVTX"));
0523     if (!tb_node)
0524     {
0525       tb_node = new PHCompositeNode("SVTX");
0526       dstNode->addNode(tb_node);
0527       if (Verbosity() > 0)
0528       {
0529         std::cout << PHWHERE << "SVTX node added" << std::endl;
0530       }
0531     }
0532 
0533     m_trackMap = new SvtxTrackMap_v1;
0534     PHIODataNode<PHObject>* tracks_node =
0535         new PHIODataNode<PHObject>(m_trackMap, m_trackMapName, "PHObject");
0536     tb_node->addNode(tracks_node);
0537     if (Verbosity() > 0)
0538     {
0539       std::cout << PHWHERE << "Svtx/" << m_trackMapName << " node added" << std::endl;
0540     }
0541   }
0542   }
0543   m_seedContainer = findNode::getClass<TrackSeedContainer>(topNode, m_trackSeedName);
0544   if (!m_seedContainer)
0545   {
0546     std::cout << PHWHERE << " Can't find track seed container " << m_trackSeedName << ", can't continue."
0547               << std::endl;
0548     return Fun4AllReturnCodes::ABORTEVENT;
0549   }
0550 
0551   m_tpcContainer = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0552   if (!m_tpcContainer)
0553   {
0554     std::cout << PHWHERE << "WARNING, TrackSeedTrackMapConverter may seg fault depending on what seeding algorithm this is run after" << std::endl;
0555   }
0556 
0557   m_siContainer = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0558   if (!m_siContainer)
0559   {
0560     std::cout << PHWHERE << "WARNING, TrackSeedTrackMapConverter may seg fault depending on what seeding algorithm this is run after" << std::endl;
0561   }
0562 
0563   m_clusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0564   if (!m_clusters)
0565   {
0566     std::cout << PHWHERE << " Can't find cluster container, can't continue."
0567               << std::endl;
0568     return Fun4AllReturnCodes::ABORTEVENT;
0569   }
0570 
0571   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0572   if (!m_tGeometry)
0573   {
0574     std::cout << PHWHERE << " Can't find ActsGeometry, can't continue."
0575               << std::endl;
0576     return Fun4AllReturnCodes::ABORTEVENT;
0577   }
0578   return Fun4AllReturnCodes::EVENT_OK;
0579 }
0580 
0581 std::pair<int, float> TrackSeedTrackMapConverter::getCosmicCharge(TrackSeed* seed,
0582                                                                   float vertexradius) const
0583 {
0584   Acts::Vector3 globalMostOuter = Acts::Vector3::Zero();
0585   Acts::Vector3 globalSecondMostOuter(0, 999999, 0);
0586   std::vector<Acts::Vector3> globpos;
0587   float largestR = 0;
0588   for (auto it = seed->begin_cluster_keys();
0589        it != seed->end_cluster_keys();
0590        ++it)
0591   {
0592     auto ckey = *it;
0593     auto clus = m_clusters->findCluster(ckey);
0594     auto glob = m_tGeometry->getGlobalPosition(ckey, clus);
0595     globpos.push_back(glob);
0596     float r = std::sqrt(square(glob.x()) + square(glob.y()));
0597     if (r > largestR && glob.y() < 0)
0598     {
0599       globalMostOuter = glob;
0600       largestR = r;
0601     }
0602   }
0603 
0604   float maxdr = std::numeric_limits<float>::max();
0605   for (auto& glob : globpos)
0606   {
0607     if (glob.y() > 0)
0608     {
0609       continue;
0610     }
0611     float dr = std::sqrt(square(globalMostOuter.x()) + square(globalMostOuter.y())) - std::sqrt(square(glob.x()) + square(glob.y()));
0612     if (dr < maxdr && dr > 10)
0613     {
0614       maxdr = dr;
0615       globalSecondMostOuter = glob;
0616     }
0617   }
0618 
0619   Acts::Vector3 vertex(0, -1 * vertexradius, 0);
0620   globalMostOuter -= vertex;
0621   globalSecondMostOuter -= vertex;
0622 
0623   const auto firstphi = atan2(globalMostOuter.y(), globalMostOuter.x());
0624   const auto secondphi = atan2(globalSecondMostOuter.y(),
0625                                globalSecondMostOuter.x());
0626   auto dphi = secondphi - firstphi;
0627   if (dphi > M_PI)
0628   {
0629     dphi = 2. * M_PI - dphi;
0630   }
0631   if (dphi < -M_PI)
0632   {
0633     dphi = 2. * M_PI + dphi;
0634   }
0635 
0636   float r1 = std::sqrt(square(globalMostOuter.x()) + square(globalMostOuter.y()));
0637   float r2 = std::sqrt(square(globalSecondMostOuter.x()) + square(globalSecondMostOuter.y()));
0638   float z1 = globalMostOuter.z();
0639   float z2 = globalSecondMostOuter.z();
0640 
0641   float slope = (r2 - r1) / (z2 - z1);
0642   int charge = 1;
0643 
0644   if (dphi > 0)
0645   {
0646     charge = -1;
0647   }
0648 
0649   return std::make_pair(charge, slope);
0650   ;
0651 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!