Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "SecondaryVertexFinder.h"
0002 #include "ActsPropagator.h"
0003 
0004 /// Tracking includes
0005 
0006 #include <trackbase/ActsGeometry.h>
0007 #include <trackbase/TpcDefs.h>
0008 #include <trackbase/TrackFitUtils.h>
0009 #include <trackbase/TrkrCluster.h>  // for TrkrCluster
0010 #include <trackbase/TrkrClusterContainer.h>
0011 #include <trackbase/TrkrDefs.h>  // for cluskey, getLayer, TrkrId
0012 
0013 #include <globalvertex/SvtxVertexMap_v1.h>
0014 #include <globalvertex/SvtxVertex_v1.h>
0015 #include <trackbase_historic/ActsTransformations.h>
0016 #include <trackbase_historic/SvtxTrack.h>  // for SvtxTrack, SvtxTrack::C...
0017 #include <trackbase_historic/SvtxTrackMap_v2.h>
0018 #include <trackbase_historic/TrackSeed_v2.h>
0019 
0020 #include <Acts/Geometry/GeometryIdentifier.hpp>
0021 #include <Acts/MagneticField/MagneticFieldProvider.hpp>
0022 #include <Acts/Surfaces/PerigeeSurface.hpp>
0023 
0024 #include <fun4all/Fun4AllReturnCodes.h>
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/PHNode.h>
0027 #include <phool/PHNodeIterator.h>
0028 
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>
0031 
0032 #include <cmath>  // for sqrt, fabs, atan2, cos
0033 #include <iomanip>
0034 #include <iostream>  // for operator<<, basic_ostream
0035 #include <map>       // for map
0036 #include <set>       // for _Rb_tree_const_iterator
0037 #include <utility>   // for pair, make_pair
0038 
0039 #include <algorithm>
0040 #include <cassert>
0041 #include <functional>
0042 #include <iostream>
0043 #include <numeric>
0044 #include <vector>
0045 
0046 #include <Eigen/Dense>
0047 
0048 #include <TFile.h>
0049 #include <TH1D.h>
0050 #include <TH2D.h>
0051 #include <TLorentzVector.h>
0052 #include <TNtuple.h>
0053 
0054 //____________________________________________________________________________..
0055 SecondaryVertexFinder::SecondaryVertexFinder(const std::string& name)
0056   : SubsysReco(name)
0057 {
0058 }
0059 
0060 //____________________________________________________________________________..
0061 SecondaryVertexFinder::~SecondaryVertexFinder() = default;
0062 
0063 //____________________________________________________________________________..
0064 int SecondaryVertexFinder::InitRun(PHCompositeNode* topNode)
0065 {
0066   if (_write_electrons_node)
0067   {
0068     CreateOutputNode(topNode);
0069   }
0070 
0071   int ret = GetNodes(topNode);
0072   if (ret != Fun4AllReturnCodes::EVENT_OK)
0073   {
0074     return ret;
0075   }
0076 
0077   if (_write_ntuple)
0078   {
0079     recomass = new TH2D("recomass", "invariant mass vs pT", 1000, 0, 5, 4000, 0, 2);
0080     hdecaypos = new TH2D("hdecaypos", "decay radius", 200, -40, 40, 200, -40, 40);
0081     hdecaypos->GetXaxis()->SetTitle("decay X (cm)");
0082     hdecaypos->GetYaxis()->SetTitle("decay Y (cm)");
0083     hdecay_radius = new TH1D("hdecay_radius", "Decay Radius", 200, 0, 40.0);
0084 
0085     ntp = new TNtuple("ntp", "decay_pairs", "x1:y1:z1:px1:py1:pz1:dca3dxy1:dca3dz1:vposx1:vposy1:vposz1:vmomx1:vmomy1:vmomz1:pca_relx_1:pca_rely_1:pca_relz_1:eta1:charge1:tpcClusters_1:quality1:eta1:x2:y2:z2:px2:py2:pz2:dca3dxy2:dca3dz2:vposx2:vposy2:vposz2:vmomx2:vmomy2:vmomz2:pca_relx_2:pca_rely_2:pca_relz_2:eta2:charge2:tpcClusters_2:quality2:eta2:vertex_x:vertex_y:vertex_z:pair_dca:invariant_mass:invariant_pt:path:has_silicon1:has_silicon2");
0086   }
0087 
0088   return ret;
0089 }
0090 
0091 void SecondaryVertexFinder::fillNtp(SvtxTrack* track1, SvtxTrack* track2, double dca3dxy1, double dca3dz1, double dca3dxy2, double dca3dz2, Eigen::Vector3d vpos1, Eigen::Vector3d vmom1, Eigen::Vector3d vpos2, Eigen::Vector3d vmom2, Acts::Vector3 pca_rel1, Acts::Vector3 pca_rel2, double pair_dca, double invariantMass, double invariantPt, double path, int has_silicon_1, int has_silicon_2)
0092 {
0093   double px1 = track1->get_px();
0094   double py1 = track1->get_py();
0095   double pz1 = track1->get_pz();
0096   auto tpcSeed1 = track1->get_tpc_seed();
0097   size_t tpcClusters1 = tpcSeed1->size_cluster_keys();
0098   double eta1 = asinh(pz1 / sqrt(pow(px1, 2) + pow(py1, 2)));
0099 
0100   double px2 = track2->get_px();
0101   double py2 = track2->get_py();
0102   double pz2 = track2->get_pz();
0103   auto tpcSeed2 = track2->get_tpc_seed();
0104   size_t tpcClusters2 = tpcSeed2->size_cluster_keys();
0105   double eta2 = asinh(pz2 / sqrt(pow(px2, 2) + pow(py2, 2)));
0106 
0107   auto vtxid = track1->get_vertex_id();
0108   auto svtxVertex = _svtx_vertex_map->get(vtxid);
0109   if (!svtxVertex)
0110   {
0111     return;
0112   }
0113 
0114   float reco_info[] = {
0115       track1->get_x(), track1->get_y(), track1->get_z(),
0116       track1->get_px(), track1->get_py(), track1->get_pz(),
0117       (float) dca3dxy1, (float) dca3dz1, (float) vpos1(0), (float) vpos1(1), (float) vpos1(2),
0118       (float) vmom1(0), (float) vmom1(1), (float) vmom1(2),
0119       (float) pca_rel1(0), (float) pca_rel1(1), (float) pca_rel1(2),
0120       (float) eta1, (float) track1->get_charge(), (float) tpcClusters1,
0121       (float) track1->get_quality(), (float) eta1,
0122       track2->get_x(), track2->get_y(), track2->get_z(),
0123       track2->get_px(), track2->get_py(), track2->get_pz(),
0124       (float) dca3dxy2, (float) dca3dz2, (float) vpos2(0), (float) vpos2(1), (float) vpos2(2),
0125       (float) vmom2(0), (float) vmom2(1), (float) vmom2(2),
0126       (float) pca_rel2(0), (float) pca_rel2(1), (float) pca_rel2(2),
0127       (float) eta2, (float) track2->get_charge(), (float) tpcClusters2,
0128       (float) track2->get_quality(), (float) eta2,
0129       svtxVertex->get_x(), svtxVertex->get_y(), svtxVertex->get_z(),
0130       (float) pair_dca, (float) invariantMass, (float) invariantPt, (float) path,
0131       (float) has_silicon_1, (float) has_silicon_2};
0132 
0133   ntp->Fill(reco_info);
0134 }
0135 
0136 //____________________________________________________________________________..
0137 int SecondaryVertexFinder::process_event(PHCompositeNode* /*topNode*/)
0138 {
0139   if (Verbosity() > 0)
0140   {
0141     std::cout << PHWHERE << " track map size " << _track_map->size() << std::endl;
0142   }
0143 
0144   // Loop over tracks and check for close DCA match with all other tracks
0145   for (auto tr1_it = _track_map->begin(); tr1_it != _track_map->end(); ++tr1_it)
0146   {
0147     auto id1 = tr1_it->first;
0148     auto tr1 = tr1_it->second;
0149 
0150     auto vertexId = tr1->get_vertex_id();
0151     const SvtxVertex* svtxVertex = _svtx_vertex_map->get(vertexId);
0152     if (!svtxVertex)
0153     {
0154       if (Verbosity() > 1)
0155       {
0156         std::cout << PHWHERE << " Failed to find vertex id " << vertexId << " skip this track " << std::endl;
0157       }
0158       continue;
0159     }
0160     if (svtxVertex->size_tracks() == 0)
0161     {
0162       continue;  // event did not have a reconstructed vertex, vertex is bogus
0163     }
0164 
0165     // Reverse or remove this to consider TPC-only tracks too
0166     if (_require_mvtx && !hasSiliconSeed(tr1))
0167     {
0168       continue;
0169     }
0170 
0171     int has_silicon_1 = 0;
0172     if (hasSiliconSeed(tr1))
0173     {
0174       has_silicon_1 = 1;
0175     }
0176 
0177     if (Verbosity() > 3)
0178     {
0179       std::cout << "Track1 " << id1 << " details: " << std::endl;
0180       outputTrackDetails(tr1);
0181     }
0182 
0183     if (tr1->get_quality() > _qual_cut)
0184     {
0185       continue;
0186     }
0187 
0188     auto tpc_seed1 = tr1->get_tpc_seed();
0189     int ntpc1 = tpc_seed1->size_cluster_keys();
0190     if (ntpc1 < 20)
0191     {
0192       continue;
0193     }
0194 
0195     float dca3dxy1, dca3dz1, dca3dxysigma1, dca3dzsigma1;
0196     get_dca(tr1, dca3dxy1, dca3dz1, dca3dxysigma1, dca3dzsigma1);
0197     // std::cout << "   tr1 " << id1 << " dca3dxy1 = " << dca3dxy1 << " dca3dz1 = " << dca3dz1 << std::endl;
0198 
0199     if (!dca3dxy1)
0200     {
0201       std::cout << " get_dca returned NAN " << std::endl;
0202       continue;
0203     }
0204     if (fabs(dca3dxy1) < _track_dcaxy_cut)
0205     {
0206       continue;
0207     }
0208     if (fabs(dca3dz1) < _track_dcaz_cut)
0209     {
0210       continue;
0211     }
0212 
0213     // look for close DCA matches with all other such tracks
0214     for (auto tr2_it = std::next(tr1_it); tr2_it != _track_map->end(); ++tr2_it)
0215     {
0216       auto id2 = tr2_it->first;
0217       auto tr2 = tr2_it->second;
0218 
0219       // Reverse or remove this to consider TPC tracks too
0220       if (_require_mvtx && !hasSiliconSeed(tr2))
0221       {
0222         continue;
0223       }
0224 
0225       int has_silicon_2 = 0;
0226       if (hasSiliconSeed(tr2))
0227       {
0228         has_silicon_2 = 1;
0229       }
0230 
0231       if (Verbosity() > 3)
0232       {
0233         std::cout << "Track2 " << id2 << " details: " << std::endl;
0234         outputTrackDetails(tr2);
0235       }
0236 
0237       if (tr2->get_charge() == tr1->get_charge())
0238       {
0239         continue;
0240       }
0241 
0242       if (tr2->get_quality() > _qual_cut)
0243       {
0244         continue;
0245       }
0246 
0247       auto tpc2_seed = tr2->get_tpc_seed();
0248       int ntpc2 = tpc2_seed->size_cluster_keys();
0249       if (ntpc2 < 20)
0250       {
0251         continue;
0252       }
0253 
0254       float dca3dxy2, dca3dz2, dca3dxysigma2, dca3dzsigma2;
0255       get_dca(tr2, dca3dxy2, dca3dz2, dca3dxysigma2, dca3dzsigma2);
0256       // std::cout << " tr2 " << id2 << " dca3dxy2 = " << dca3dxy2 << " dca3dz2 = " << dca3dz2 << std::endl;
0257       if (!dca3dxy2)
0258       {
0259         std::cout << " get_dca returned NAN " << std::endl;
0260         continue;
0261       }
0262       if (fabs(dca3dxy2) < _track_dcaxy_cut)
0263       {
0264         continue;
0265       }
0266       if (fabs(dca3dz2) < _track_dcaz_cut)
0267       {
0268         continue;
0269       }
0270 
0271       // find DCA and PCA of these two tracks
0272       if (Verbosity() > 3)
0273       {
0274         std::cout << "Check pair DCA for tracks " << id1 << " and  " << id2 << std::endl;
0275       }
0276 
0277       double R1;
0278       Eigen::Vector2d center1;
0279       getCircleXYTrack(tr1, R1, center1);
0280       if (Verbosity() > 2)
0281       {
0282         std::cout << std::endl
0283                   << "Track 1 circle: " << std::endl;
0284         std::cout << "  tr1.x " << tr1->get_x() << " tr1.y " << tr1->get_y() << " tr1.z " << tr1->get_z() << " charge " << tr1->get_charge() << std::endl;
0285         std::cout << "   tr1.px " << tr1->get_px() << " tr1.py " << tr1->get_py() << " tr1.pz " << tr1->get_pz() << std::endl;
0286         std::cout << "         track1 " << tr1->get_id() << " circle R " << R1 << " (x, y)  " << center1(0) << "  " << center1(1) << std::endl;
0287       }
0288       double R2;
0289       Eigen::Vector2d center2;
0290       getCircleXYTrack(tr2, R2, center2);
0291       if (Verbosity() > 2)
0292       {
0293         std::cout << "Track 2 circle: " << std::endl;
0294         std::cout << "   tr2.x " << tr2->get_x() << " tr2.y " << tr2->get_y() << " tr2.z " << tr2->get_z() << " charge " << tr2->get_charge() << std::endl;
0295         std::cout << "   tr2.px " << tr2->get_px() << " tr2.py " << tr2->get_py() << " tr2.pz " << tr2->get_pz() << std::endl;
0296         std::cout << "         track2 " << tr2->get_id() << " circle R " << R2 << " (x, y)  " << center2(0) << "  " << center2(1) << std::endl;
0297       }
0298 
0299       std::vector<double> intersections;
0300       if (!circle_circle_intersection(R1, center1(0), center1(1), R2, center2(0), center2(1), intersections))
0301       {
0302         if (Verbosity() > 2)
0303         {
0304           std::cout << "    - no intersections, skip this pair" << std::endl;
0305         }
0306         continue;
0307       }
0308 
0309       Eigen::Vector2d intersection[2];
0310       if (intersections.size() == 2)
0311       {
0312         intersection[0](0) = intersections[0];
0313         intersection[0](1) = intersections[1];
0314       }
0315       if (intersections.size() == 4)
0316       {
0317         intersection[1](0) = intersections[2];
0318         intersection[1](1) = intersections[3];
0319       }
0320 
0321       // process both intersections
0322       for (int i = 0; i < 2; ++i)
0323       {
0324         if (intersection[i].norm() == 0)
0325         {
0326           continue;
0327         }
0328 
0329         double vradius = sqrt(intersection[i](0) * intersection[i](0) + intersection[i](1) * intersection[i](1));
0330         if (vradius > _max_intersection_radius)
0331         {
0332           continue;
0333         }
0334 
0335         double z1 = getZFromIntersectionXY(tr1, R1, center1, intersection[i]);
0336         double z2 = getZFromIntersectionXY(tr2, R2, center2, intersection[i]);
0337         if (Verbosity() > 2)
0338         {
0339           std::cout << " track intersection " << i << " at (x,y) " << intersection[i](0) << "  " << intersection[i](1)
0340                     << " radius " << vradius << " est. z1 " << z1 << " est. z2 " << z2 << std::endl;
0341         }
0342 
0343         if (fabs(z1 - z2) > 2.0)  // skip this intersection, it is not good
0344         {
0345           if (Verbosity() > 2)
0346           {
0347             std::cout << "    z-mismatch - wrong intersection, skip it " << std::endl;
0348           }
0349           continue;
0350         }
0351 
0352         Eigen::Vector3d vpos1(0, 0, 0), vmom1(0, 0, 0);
0353         Eigen::Vector3d vpos2(0, 0, 0), vmom2(0, 0, 0);
0354 
0355         // Project the tracks to the intersection
0356         Eigen::Vector3d intersect1(intersection[i](0), intersection[i](1), z1);
0357         if (!projectTrackToPoint(tr1, intersect1, vpos1, vmom1))
0358         {
0359           continue;
0360         }
0361         if (Verbosity() > 2)
0362         {
0363           std::cout << "  Projected track 1 to point " << intersect1(0) << "  " << intersect1(1) << "  " << intersect1(2) << std::endl;
0364           std::cout << "                                 has vpos " << vpos1(0) << "  " << vpos1(1) << "  " << vpos1(2) << std::endl;
0365         }
0366         Eigen::Vector3d intersect2(intersection[i](0), intersection[i](1), z2);
0367         if (!projectTrackToPoint(tr2, intersect2, vpos2, vmom2))
0368         {
0369           continue;
0370         }
0371         if (Verbosity() > 2)
0372         {
0373           std::cout << "  Projected track 2 to point " << intersect2(0) << "  " << intersect2(1) << "  " << intersect2(2) << std::endl;
0374           std::cout << "                                 has vpos " << vpos2(0) << "  " << vpos2(1) << "  " << vpos2(2) << std::endl;
0375         }
0376 
0377         // check that the z positions are close
0378         if (fabs(vpos1(2) - vpos2(2)) > _projected_track_z_cut)
0379         {
0380           if (Verbosity() > 0)
0381           {
0382             std::cout << "    Warning: projected z positions are screwed up, should not be" << std::endl;
0383           }
0384           continue;
0385         }
0386 
0387         if (Verbosity() > 2)
0388         {
0389           std::cout << "Summary for projected pair:" << std::endl;
0390           std::cout << " Fitted tracks: " << std::endl;
0391           std::cout << "   tr1.x " << tr1->get_x() << " tr1.y " << tr1->get_y() << " tr1.z " << tr1->get_z() << std::endl;
0392           std::cout << "   tr2.x " << tr2->get_x() << " tr2.y " << tr2->get_y() << " tr2.z " << tr2->get_z() << std::endl;
0393           std::cout << "   tr1.px " << tr1->get_px() << " tr1.py " << tr1->get_py() << " tr1.pz " << tr1->get_pz() << std::endl;
0394           std::cout << "   tr2.px " << tr2->get_px() << " tr2.py " << tr2->get_py() << " tr2.pz " << tr2->get_pz() << std::endl;
0395           std::cout << " Projected tracks: " << std::endl;
0396           std::cout << "   pos1.x " << vpos1(0) << " pos1.y " << vpos1(1) << " pos1.z " << vpos1(2) << std::endl;
0397           std::cout << "   pos2.x " << vpos2(0) << " pos2.y " << vpos2(1) << " pos2.z " << vpos2(2) << std::endl;
0398           std::cout << "   mom1.x " << vmom1(0) << " mom1.y " << vmom1(1) << " mom1.z " << vmom1(2) << std::endl;
0399           std::cout << "   mom2.x " << vmom2(0) << " mom2.y " << vmom2(1) << " mom2.z " << vmom2(2) << std::endl;
0400         }
0401 
0402         // Improve the pair dca using a local straight line approximation
0403         double pair_dca;
0404         Eigen::Vector3d PCA1(0, 0, 0), PCA2(0, 0, 0);
0405         findPcaTwoLines(vpos1, vmom1, vpos2, vmom2, pair_dca, PCA1, PCA2);
0406         if (Verbosity() > 2)
0407         {
0408           std::cout << "       pair_dca " << pair_dca << " two_track_dcacut " << _two_track_dcacut << std::endl;
0409           std::cout << "       PCA1 " << PCA1(0) << "  " << PCA1(1) << "  " << PCA1(2) << std::endl;
0410           std::cout << "       PCA2 " << PCA2(0) << "  " << PCA2(1) << "  " << PCA2(2) << std::endl;
0411         }
0412 
0413         if (fabs(pair_dca) > _two_track_dcacut)
0414         {
0415           continue;
0416         }
0417 
0418         // calculate the invariant mass using the track states at the decay vertex
0419 
0420         TLorentzVector t1;
0421         Float_t E1 = sqrt(pow(vmom1(0), 2) + pow(vmom1(1), 2) + pow(vmom1(2), 2) + pow(_decaymass, 2));
0422         t1.SetPxPyPzE(vmom1(0), vmom1(1), vmom1(2), E1);
0423 
0424         TLorentzVector t2;
0425         Float_t E2 = sqrt(pow(vmom2(0), 2) + pow(vmom2(1), 2) + pow(vmom2(2), 2) + pow(_decaymass, 2));
0426         t2.SetPxPyPzE(vmom2(0), vmom2(1), vmom2(2), E2);
0427 
0428         TLorentzVector tsum = t1 + t2;
0429 
0430         // calculate the decay length
0431         Eigen::Vector3d PCA = (vpos1 + vpos2) / 2.0;  // average the PCA of the track pair
0432         auto vtxid = tr1->get_vertex_id();
0433         auto vertex1 = _svtx_vertex_map->get(vtxid);
0434         Eigen::Vector3d VTX(vertex1->get_x(), vertex1->get_y(), vertex1->get_z());
0435         Eigen::Vector3d path = PCA - VTX;
0436         double decay_radius = sqrt(pow(PCA(0), 2) + pow(PCA(1), 2));
0437 
0438         if (path.norm() > _min_path_cut)
0439         {
0440           if (Verbosity() > 0)
0441           {
0442             std::cout << "    Pair mass " << tsum.M() << " pair pT " << tsum.Pt()
0443                       << " decay length " << path.norm() << std::endl;
0444           }
0445 
0446           if (_write_ntuple)
0447           {
0448             fillNtp(tr1, tr2, dca3dxy1, dca3dz1, dca3dxy2, dca3dz2, vpos1, vmom1, vpos2, vmom2,
0449                     PCA1, PCA2, pair_dca, tsum.M(), tsum.Pt(), path.norm(), has_silicon_1, has_silicon_2);
0450           }
0451 
0452           if (_write_electrons_node)
0453           {
0454             if (passConversionElectronCuts(tsum, tr1, tr2, pair_dca, PCA, VTX))
0455             {
0456               if (Verbosity() > 0)
0457               {
0458                 std::cout << "     **** inserting tracks " << tr1->get_id() << "  and " << tr2->get_id() << std::endl;
0459               }
0460 
0461               // Add decay particles to output node
0462               _track_map_electrons->insertWithKey(tr1, tr1->get_id());
0463               _track_map_electrons->insertWithKey(tr2, tr2->get_id());
0464 
0465               if (_write_ntuple)
0466               {
0467                 // these are just to check on the effect of the cuts
0468                 recomass->Fill(tsum.Pt(), tsum.M());
0469                 hdecay_radius->Fill(decay_radius);
0470                 hdecaypos->Fill(PCA(0), PCA(1));
0471               }
0472             }
0473           }
0474         }
0475       }
0476     }
0477   }
0478 
0479   if (Verbosity() > 0)
0480   {
0481     std::cout << PHWHERE << " electron track map size " << _track_map_electrons->size() << std::endl;
0482     for (auto& _track_map_electron : *_track_map_electrons)
0483     {
0484       auto id = _track_map_electron.first;
0485       auto tr = _track_map_electron.second;
0486 
0487       std::cout << " Electron track " << id
0488                 << " x " << tr->get_x() << " y " << tr->get_y() << " z " << tr->get_z()
0489                 << " px " << tr->get_px() << " py " << tr->get_py() << " pz " << tr->get_pz()
0490                 << std::endl;
0491     }
0492   }
0493 
0494   return Fun4AllReturnCodes::EVENT_OK;
0495 }
0496 
0497 bool SecondaryVertexFinder::passConversionElectronCuts(const TLorentzVector& tsum, SvtxTrack* tr1, SvtxTrack* tr2, float pair_dca, const Eigen::Vector3d& PCA, const Eigen::Vector3d& VTX)
0498 {
0499   bool pass = false;
0500 
0501   // additional single track cuts
0502   auto tpc_seed1 = tr1->get_tpc_seed();
0503   unsigned int ntpc1 = tpc_seed1->size_cluster_keys();
0504   auto tpc_seed2 = tr2->get_tpc_seed();
0505   unsigned int ntpc2 = tpc_seed2->size_cluster_keys();
0506   if (ntpc1 < _min_tpc_clusters || ntpc2 < _min_tpc_clusters)
0507   {
0508     return pass;
0509   }
0510 
0511   // pair cuts
0512   if (fabs(pair_dca) > _conversion_pair_dcacut)
0513   {
0514     return pass;
0515   };
0516 
0517   if (tsum.M() > _max_mass_cut)
0518   {
0519     return pass;
0520   }
0521 
0522   if (tsum.Pt() < _invariant_pt_cut)
0523   {
0524     return pass;
0525   }
0526 
0527   Eigen::Vector3d path = PCA - VTX;
0528   if (path.norm() < 0.2)
0529   {
0530     return pass;
0531   }
0532 
0533   // Angle between path vector and reco pair momentum vector
0534   Eigen::Vector3d invariant_mom(tsum.Px(), tsum.Py(), tsum.Pz());
0535   double costheta = path.dot(invariant_mom) / (path.norm() * invariant_mom.norm());
0536   //      std::cout << " costheta " << costheta << std::endl;
0537   if (costheta < _costheta_cut)
0538   {
0539     return pass;
0540   };
0541 
0542   // selects "decays" with zero mass
0543   if (fabs(tr1->get_eta() - tr2->get_eta()) > _deta_cut)
0544   {
0545     return pass;
0546   };
0547 
0548   pass = true;
0549   return pass;
0550 }
0551 
0552 double SecondaryVertexFinder::getZFromIntersectionXY(SvtxTrack* track, double& R, Eigen::Vector2d& center, Eigen::Vector2d intersection)
0553 {
0554   // Starting at the track base position, the path in the XY plane and the path in the Z direction are proportional
0555   // They are related by:   dpath/dz = pT/pz
0556 
0557   // The xy path to the intersection is an arc due to a rotation of pT:
0558   double phi0 = atan2(track->get_y() - center(1), track->get_x() - center(0));
0559   double phi1 = atan2(intersection(1) - center(1), intersection(0) - center(0));
0560   double xypath = R * fabs(phi1 - phi0);
0561 
0562   double zpath = xypath * track->get_pz() / track->get_pt();
0563   double z = track->get_z() + zpath;
0564   if (Verbosity() > 3)
0565   {
0566     std::cout << " Radius " << R << " center xy " << center(0) << "  " << center(1) << "  pT " << track->get_pt() << " pz " << track->get_pz() << std::endl;
0567     std::cout << " phi0 " << phi0 << " y0 " << track->get_y() << " x0 " << track->get_x() << std::endl;
0568     std::cout << " phi1 " << phi1 << " intersection y " << intersection(1) << " intersection x " << intersection(0) << std::endl;
0569     std::cout << " xypath " << xypath << " zpath " << zpath << " z0 " << track->get_z() << " z " << z << std::endl;
0570   }
0571 
0572   return z;
0573 }
0574 
0575 void SecondaryVertexFinder::getCircleXYTrack(SvtxTrack* track, double& R, Eigen::Vector2d& center)
0576 {
0577   // Get the circle equation for the track
0578   double Bz = 1.4;  // T
0579   double x = track->get_x();
0580   double y = track->get_y();
0581   double px = track->get_px();
0582   double py = track->get_py();
0583   double pt = track->get_pt();
0584   int charge = track->get_charge();
0585 
0586   // radius of curvature is from pT and B field
0587   R = 3.3 * pt / (Bz * (double) charge) * 100;  // convert from m to cm, sign changes for negative charge
0588 
0589   // the normal unit vector in (x,y) space at (x,y) is:
0590   Eigen::Vector2d normal(py, -px);
0591   normal /= normal.norm();
0592 
0593   // The circle center is at
0594   Eigen::Vector2d base(x, y);
0595   center = base + R * normal;
0596 
0597   // needed the sign to get the direction of the center, now we drop it
0598   R = fabs(R);
0599 
0600   return;
0601 }
0602 
0603 bool SecondaryVertexFinder::projectTrackToPoint(SvtxTrack* track, Eigen::Vector3d& PCA, Eigen::Vector3d& pos, Eigen::Vector3d& mom)
0604 {
0605   bool ret = true;
0606   PCA *= Acts::UnitConstants::cm;
0607 
0608   /// create perigee surface
0609   auto perigee = Acts::Surface::makeShared<Acts::PerigeeSurface>(PCA);
0610 
0611   ActsPropagator propagator(_tGeometry);
0612   auto params = propagator.makeTrackParams(track, _svtx_vertex_map);
0613   if (!params.ok())
0614   {
0615     /// Couldn't construct strict PCA track parameter requirement
0616     return false;
0617   }
0618 
0619   propagator.verbosity(Verbosity());
0620   auto result = propagator.propagateTrack(params.value(), perigee);
0621   if (result.ok())
0622   {
0623     auto resultparams = result.value().second;
0624     auto projectionPos = resultparams.position(_tGeometry->geometry().getGeoContext());
0625     const auto momentum = resultparams.momentum();
0626     pos(0) = projectionPos.x() / Acts::UnitConstants::cm;
0627     pos(1) = projectionPos.y() / Acts::UnitConstants::cm;
0628     pos(2) = projectionPos.z() / Acts::UnitConstants::cm;
0629 
0630     mom(0) = momentum.x();
0631     mom(1) = momentum.y();
0632     mom(2) = momentum.z();
0633   }
0634   else
0635   {
0636     if (Verbosity() > 0)
0637     {
0638       std::cout << "   Failed projection of track with: " << std::endl;
0639       std::cout << " x,y,z = " << track->get_x() << "  " << track->get_y() << "  " << track->get_z() << std::endl;
0640       std::cout << " px,py,pz = " << track->get_px() << "  " << track->get_py() << "  " << track->get_pz() << std::endl;
0641       std::cout << " to point (x,y,z) = " << PCA(0) / Acts::UnitConstants::cm << "  "
0642                 << PCA(1) / Acts::UnitConstants::cm << "  " << PCA(2) / Acts::UnitConstants::cm << std::endl;
0643     }
0644 
0645     ret = false;
0646   }
0647 
0648   return ret;
0649 }
0650 
0651 void SecondaryVertexFinder::outputTrackDetails(SvtxTrack* tr)
0652 {
0653   auto tpc_seed = tr->get_tpc_seed();
0654   int ntpc = tpc_seed->size_cluster_keys();
0655 
0656   auto silicon_seed = tr->get_silicon_seed();
0657 
0658   int nsilicon = 0;
0659   if (silicon_seed)
0660   {
0661     nsilicon = silicon_seed->size_cluster_keys();
0662   }
0663 
0664   auto pt = tr->get_pt();
0665   auto eta = tr->get_eta();
0666   auto x = tr->get_x();
0667   auto y = tr->get_y();
0668   auto z = tr->get_z();
0669   auto qual = tr->get_quality();
0670 
0671   std::cout << "   ntpc " << ntpc << " nsilicon " << nsilicon << " quality " << qual
0672             << " eta " << eta << std::endl;
0673   std::cout << "   pt " << pt << " x " << x << " y " << y << " z " << z << std::endl;
0674 
0675   auto vtxid = tr->get_vertex_id();
0676   auto vertex = _svtx_vertex_map->get(vtxid);
0677   std::cout << "    vtxid " << vtxid
0678             << " vertex x " << vertex->get_x()
0679             << " vertex y " << vertex->get_y()
0680             << " vertex z " << vertex->get_z()
0681             << std::endl;
0682 }
0683 
0684 bool SecondaryVertexFinder::hasSiliconSeed(SvtxTrack* tr)
0685 {
0686   bool ret = false;
0687   auto silicon_seed = tr->get_silicon_seed();
0688   if (silicon_seed)
0689   {
0690     ret = true;
0691   }
0692 
0693   return ret;
0694 }
0695 
0696 void SecondaryVertexFinder::get_dca(SvtxTrack* track,
0697                                     float& dca3dxy, float& dca3dz,
0698                                     float& dca3dxysigma, float& dca3dzsigma)
0699 {
0700   dca3dxy = NAN;
0701   Acts::Vector3 pos(track->get_x(),
0702                     track->get_y(),
0703                     track->get_z());
0704   Acts::Vector3 mom(track->get_px(),
0705                     track->get_py(),
0706                     track->get_pz());
0707 
0708   auto vtxid = track->get_vertex_id();
0709   auto svtxVertex = _svtx_vertex_map->get(vtxid);
0710   if (!svtxVertex)
0711   {
0712     std::cout << "   Failed to find vertex for track " << std::endl;
0713     return;
0714   }
0715 
0716   Acts::Vector3 vertex(svtxVertex->get_x(),
0717                        svtxVertex->get_y(),
0718                        svtxVertex->get_z());
0719 
0720   if (Verbosity() > 3)
0721   {
0722     std::cout << "   track " << track->get_id() << " vertex id is " << vtxid
0723               << " vertex is " << svtxVertex->get_x() << "  "
0724               << svtxVertex->get_y() << "  "
0725               << svtxVertex->get_z() << std::endl;
0726   }
0727 
0728   pos -= vertex;
0729 
0730   Acts::ActsSquareMatrix<3> posCov;
0731   for (int i = 0; i < 3; ++i)
0732   {
0733     for (int j = 0; j < 3; ++j)
0734     {
0735       posCov(i, j) = track->get_error(i, j);
0736     }
0737   }
0738 
0739   Acts::Vector3 r = mom.cross(Acts::Vector3(0., 0., 1.));
0740   float phi = atan2(r(1), r(0));
0741   phi *= -1.0;  // rotates vector clockwise to x axis
0742 
0743   Acts::RotationMatrix3 rot;
0744   Acts::RotationMatrix3 rot_T;
0745   rot(0, 0) = cos(phi);
0746   rot(0, 1) = -sin(phi);
0747   rot(0, 2) = 0;
0748   rot(1, 0) = sin(phi);
0749   rot(1, 1) = cos(phi);
0750   rot(1, 2) = 0;
0751   rot(2, 0) = 0;
0752   rot(2, 1) = 0;
0753   rot(2, 2) = 1;
0754 
0755   rot_T = rot.transpose();
0756 
0757   Acts::Vector3 pos_R = rot * pos;
0758   Acts::ActsSquareMatrix<3> rotCov = rot * posCov * rot_T;
0759 
0760   dca3dxy = pos_R(0);
0761   dca3dz = pos_R(2);
0762   dca3dxysigma = sqrt(rotCov(0, 0));
0763   dca3dzsigma = sqrt(rotCov(2, 2));
0764 }
0765 
0766 void SecondaryVertexFinder::findPcaTwoLines(Eigen::Vector3d pos1, Eigen::Vector3d mom1, Eigen::Vector3d pos2, Eigen::Vector3d mom2,
0767                                             double& dca, Eigen::Vector3d& PCA1, Eigen::Vector3d& PCA2)
0768 {
0769   Eigen::Vector3d a1(pos1(0), pos1(1), pos1(2));
0770   Eigen::Vector3d b1(mom1(0) / mom1.norm(), mom1(1) / mom1.norm(), mom1(2) / mom1.norm());
0771   Eigen::Vector3d a2(pos2(0), pos2(1), pos2(2));
0772   Eigen::Vector3d b2(mom2(0) / mom2.norm(), mom2(1) / mom2.norm(), mom2(2) / mom2.norm());
0773 
0774   // The shortest distance between two skew lines described by
0775   //  a1 + c * b1,   a2 + d * b2
0776   // where a1, a2, are vectors representing points on the lines, b1, b2 are direction vectors,
0777   //  and c and d are scalars, is:
0778   // dca = (b1 x b2) .(a2-a1) / |b1 x b2|
0779 
0780   // bcrossb/mag_bcrossb is a unit vector perpendicular to both direction vectors b1 and b2
0781   auto bcrossb = b1.cross(b2);
0782   auto mag_bcrossb = bcrossb.norm();
0783   // a2-a1 is the vector joining any arbitrary points on the two lines
0784   auto aminusa = a2 - a1;
0785 
0786   // The DCA of these two lines is the projection of a2-a1 along the direction of the perpendicular to both
0787   // remember that a2-a1 is longer than (or equal to) the dca by definition
0788   dca = 999;
0789   if (mag_bcrossb != 0)
0790   {
0791     dca = bcrossb.dot(aminusa) / mag_bcrossb;
0792   }
0793   else
0794   {
0795     return;  // same track, skip combination
0796   }
0797 
0798   // get the points at which the normal to the lines intersect the lines, where the lines are perpendicular
0799 
0800   double X = b1.dot(b2) - b1.dot(b1) * b2.dot(b2) / b2.dot(b1);
0801   double Y = (a2.dot(b2) - a1.dot(b2)) - (a2.dot(b1) - a1.dot(b1)) * b2.dot(b2) / b2.dot(b1);
0802   double c = Y / X;
0803 
0804   double F = b1.dot(b1) / b2.dot(b1);
0805   double G = -(a2.dot(b1) - a1.dot(b1)) / b2.dot(b1);
0806   double d = c * F + G;
0807 
0808   // then the points of closest approach are:
0809   PCA1 = a1 + c * b1;
0810   PCA2 = a2 + d * b2;
0811 
0812   return;
0813 }
0814 
0815 //===========================
0816 // replace with a call to TrackFitUtils
0817 //===========================
0818 bool SecondaryVertexFinder::circle_circle_intersection(double r0, double x0, double y0, double r1, double x1, double y1, std::vector<double>& intersectionXY)
0819 {
0820   bool ret = true;
0821 
0822   Eigen::Vector2d p0(x0, y0);
0823   Eigen::Vector2d p1(x1, y1);
0824 
0825   double d = (p0 - p1).norm();
0826 
0827   //  std::cout << "   d " << d << " r1-r0 " << fabs(r1-r0) << std::endl;
0828 
0829   // no intersection possible
0830   if (d < fabs(r1 - r0))
0831   {
0832     return false;  // one circle inside the other
0833   }
0834   if (d > r0 + r1)
0835   {
0836     // careful: conversion electrons will look like zero mass decays
0837     // fluctuations may cause the circles to (just) not touch - what to do about that?
0838     //   if d - (r0+r1) < dr,   then there is only one PCA, and it is on the line between the two centers
0839     double dr = 0.2;  // 2 mm
0840     if (fabs(d - (r0 + r1)) < dr)
0841     {
0842       // find the closest point on circle 0 to the center of circle 1
0843       Eigen::Vector2d u0 = (p1 - p0);
0844       u0 /= u0.norm();
0845       Eigen::Vector2d PCA0 = p0 + u0 * r0;
0846 
0847       Eigen::Vector2d u1 = (p0 - p1);
0848       u1 /= u1.norm();
0849       Eigen::Vector2d PCA1 = p1 + u1 * r1;
0850 
0851       auto PCA = (PCA0 + PCA1) / 2.0;
0852       intersectionXY.push_back(PCA(0));
0853       intersectionXY.push_back(PCA(1));
0854 
0855       if (Verbosity() > 2)
0856       {
0857         std::cout << "      *** Special case: Barely touching circles: "
0858                   << " PCA.x, PCA.y " << PCA(0) << "   " << PCA(1) << std::endl;
0859       }
0860       return ret;
0861     }
0862     else
0863     {
0864       return false;
0865     }
0866   }
0867 
0868   double a = (r0 * r0 - r1 * r1 + d * d) / (2 * d);
0869   double h = sqrt(r0 * r0 - a * a);
0870 
0871   double x2 = x0 + a * (x1 - x0) / d;
0872   double y2 = y0 + a * (y1 - y0) / d;
0873 
0874   double x3a = x2 + h * (y1 - y0) / d;  // also x3=x2-h*(y1-y0)/d
0875   double y3a = y2 - h * (x1 - x0) / d;  // also y3=y2+h*(x1-x0)/d
0876 
0877   double x3b = x2 - h * (y1 - y0) / d;  // also x3=x2-h*(y1-y0)/d
0878   double y3b = y2 + h * (x1 - x0) / d;  // also y3=y2+h*(x1-x0)/d
0879 
0880   intersectionXY.push_back(x3a);
0881   intersectionXY.push_back(y3a);
0882   intersectionXY.push_back(x3b);
0883   intersectionXY.push_back(y3b);
0884 
0885   return ret;
0886 }
0887 
0888 int SecondaryVertexFinder::End(PHCompositeNode* /*topNode*/)
0889 {
0890   if (_write_ntuple)
0891   {
0892     TFile* fout = new TFile(outfile.c_str(), "recreate");
0893     recomass->Write();
0894     hdecaypos->Write();
0895     hdecay_radius->Write();
0896     ntp->Write();
0897     fout->Close();
0898   }
0899 
0900   return Fun4AllReturnCodes::EVENT_OK;
0901 }
0902 
0903 int SecondaryVertexFinder::CreateOutputNode(PHCompositeNode* topNode)
0904 {
0905   PHNodeIterator iter(topNode);
0906 
0907   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0908 
0909   if (!dstNode)
0910   {
0911     std::cerr << "DST node is missing, quitting" << std::endl;
0912     throw std::runtime_error("Failed to find DST node in SecondaryVertexFinder");
0913   }
0914 
0915   PHNodeIterator dstIter(topNode);
0916   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0917 
0918   if (!svtxNode)
0919   {
0920     svtxNode = new PHCompositeNode("SVTX");
0921     dstNode->addNode(svtxNode);
0922   }
0923 
0924   if (_write_electrons_node)
0925   {
0926     _track_map_electrons = new SvtxTrackMap_v2;
0927     auto tracks_node = new PHIODataNode<PHObject>(_track_map_electrons, "SvtxTrackMapElectrons", "PHObject");
0928     svtxNode->addNode(tracks_node);
0929     if (Verbosity())
0930     {
0931       std::cout << "Svtx/SvtxTrackMapElectrons node added" << std::endl;
0932     }
0933   }
0934 
0935   return true;
0936 }
0937 
0938 int SecondaryVertexFinder::GetNodes(PHCompositeNode* topNode)
0939 {
0940   _track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0941   if (!_track_map)
0942   {
0943     std::cout << PHWHERE << " ERROR: Can't find SvtxTrackMap: " << std::endl;
0944     return Fun4AllReturnCodes::ABORTEVENT;
0945   }
0946 
0947   /*
0948   _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0949   if (!_cluster_map)
0950     {
0951       std::cout << PHWHERE << " ERROR: Can't find node TRKR_CLUSTER" << std::endl;
0952       return Fun4AllReturnCodes::ABORTEVENT;
0953     }
0954   */
0955 
0956   _svtx_vertex_map = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0957   if (!_svtx_vertex_map)
0958   {
0959     std::cout << PHWHERE << "No vertex node, quit!" << std::endl;
0960     return Fun4AllReturnCodes::ABORTEVENT;
0961   }
0962 
0963   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0964   if (!_tGeometry)
0965   {
0966     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0967     return Fun4AllReturnCodes::ABORTEVENT;
0968   }
0969 
0970   return Fun4AllReturnCodes::EVENT_OK;
0971 }
[ 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!