Back to home page

sPhenix code displayed by LXR
 
 

    

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

0001 #include "KshortReconstruction.h"
0002 
0003 #include <ffaobjects/EventHeader.h>
0004 
0005 #include <trackbase/TrkrDefs.h>
0006 #include <trackbase/ActsGeometry.h>
0007 
0008 #include <trackbase_historic/ActsTransformations.h>
0009 #include <trackbase_historic/SvtxTrackMap_v2.h>
0010 
0011 #include <trackreco/ActsPropagator.h>
0012 
0013 #include <globalvertex/SvtxVertexMap.h>
0014 
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 
0017 #include <phool/PHCompositeNode.h>
0018 #include <phool/getClass.h>
0019 
0020 #include <Acts/Surfaces/CylinderSurface.hpp>
0021 
0022 #pragma GCC diagnostic push
0023 #pragma GCC diagnostic ignored "-Wshadow"
0024 #include <ActsExamples/EventData/Trajectories.hpp>
0025 #pragma GCC diagnostic pop
0026 
0027 #include <TFile.h>
0028 #include <TH1.h>
0029 #include <TLorentzVector.h>
0030 #include <TNtuple.h>
0031 #include <TSystem.h>
0032 
0033 #include <cmath>
0034 #include <utility>
0035 
0036 using BoundTrackParam = const Acts::BoundTrackParameters;
0037 using BoundTrackParamResult = Acts::Result<BoundTrackParam>;
0038 using SurfacePtr = std::shared_ptr<const Acts::Surface>;
0039 using Trajectory = ActsExamples::Trajectories;
0040 
0041 
0042 int KshortReconstruction::process_event(PHCompositeNode* topNode)
0043 {
0044 
0045   EventHeader* evtHeader = findNode::getClass<EventHeader>(topNode, "EventHeader");
0046 
0047   int m_runNumber;
0048   int m_evtNumber;
0049   if (evtHeader)
0050   {
0051     m_runNumber = evtHeader->get_RunNumber();
0052     m_evtNumber = evtHeader->get_EvtSequence();
0053   }
0054   else
0055   {
0056     m_runNumber = m_evtNumber = -1;
0057   }
0058 
0059   // Loop over tracks and check for close DCA match with all other tracks
0060   for (auto tr1_it = m_svtxTrackMap->begin(); tr1_it != m_svtxTrackMap->end(); ++tr1_it)
0061   {
0062     auto id1 = tr1_it->first;
0063     auto *tr1 = tr1_it->second;
0064     if (tr1->get_quality() > _qual_cut)
0065     {
0066       continue;
0067     }
0068     if (tr1->get_pt() < track_pt_cut)
0069     {
0070       continue;
0071     }
0072 
0073     short int crossing1 = tr1->get_crossing();
0074 
0075     // calculate number silicon tracks
0076     double this_dca_cut = track_dca_cut;
0077     TrackSeed* siliconseed = tr1->get_silicon_seed();
0078 
0079     if (!siliconseed)
0080     {
0081       this_dca_cut *= 5;
0082       if (Verbosity() > 2)
0083       {
0084         std::cout << "silicon seed not found" << std::endl;
0085       }
0086       if (_require_mvtx)
0087       {
0088         continue;
0089       }
0090     }
0091 
0092     std::vector<unsigned int> nstates1 = getTrackStates(tr1);
0093     unsigned int track1_mvtx_state_size = nstates1[0];
0094     unsigned int track1_intt_state_size = nstates1[1];
0095     // unsigned int track1_tpc_state_size = nstates1[2];
0096     // unsigned int track1_mms_state_size = nstates1[3];
0097     
0098     unsigned int track1_silicon_cluster_size = std::numeric_limits<unsigned int>::quiet_NaN();
0099     if (siliconseed)
0100     {
0101       track1_silicon_cluster_size = siliconseed->size_cluster_keys();
0102     }
0103 
0104     std::vector<TrkrDefs::cluskey> ckeys1;
0105     if (siliconseed)
0106     {
0107       ckeys1.insert(ckeys1.end(), siliconseed->begin_cluster_keys(), siliconseed->end_cluster_keys());
0108     }
0109 
0110     unsigned int track1_mvtx_cluster_size = 0;
0111     unsigned int track1_intt_cluster_size = 0;
0112     for (const auto& ckey : ckeys1)
0113     {
0114       auto detid = TrkrDefs::getTrkrId(ckey);
0115       if (detid == TrkrDefs::TrkrId::mvtxId)
0116       {
0117         track1_mvtx_cluster_size++;
0118       }
0119       else if (detid == TrkrDefs::TrkrId::inttId)
0120       {
0121         track1_intt_cluster_size++;
0122       }
0123     }
0124 
0125     Acts::Vector3 pos1(tr1->get_x(), tr1->get_y(), tr1->get_z());
0126     Acts::Vector3 mom1(tr1->get_px(), tr1->get_py(), tr1->get_pz());
0127     Acts::Vector3 dcaVals1 = calculateDca(tr1, mom1, pos1);
0128     // first dca cuts
0129     if (fabs(dcaVals1(0)) < this_dca_cut or fabs(dcaVals1(1)) < this_dca_cut)
0130     {
0131       continue;
0132     }
0133 
0134     // look for close DCA matches with all other such tracks
0135     for (auto tr2_it = std::next(tr1_it); tr2_it != m_svtxTrackMap->end(); ++tr2_it)
0136     {
0137       auto id2 = tr2_it->first;
0138       auto *tr2 = tr2_it->second;
0139       if (tr2->get_quality() > _qual_cut)
0140       {
0141         continue;
0142       }
0143       if (tr2->get_pt() < track_pt_cut)
0144       {
0145         continue;
0146       }
0147 
0148       short int crossing2 = tr2->get_crossing();
0149 
0150       // calculate number silicon tracks
0151       TrackSeed* siliconseed2 = tr2->get_silicon_seed();
0152 
0153       double this_dca_cut2 = track_dca_cut;
0154 
0155       if (!siliconseed2)
0156       {
0157         this_dca_cut2 *= 5;
0158         if (Verbosity() > 2)
0159         {
0160           std::cout << "silicon seed not found" << std::endl;
0161         }
0162         if (_require_mvtx)
0163         {
0164           continue;
0165         }
0166       }
0167 
0168       std::vector<unsigned int> nstates2 = getTrackStates(tr2);
0169       unsigned int track2_mvtx_state_size = nstates2[0];
0170       unsigned int track2_intt_state_size = nstates2[1];
0171       // unsigned int track2_tpc_state_size = nstates2[2];
0172       // unsigned int track2_mms_state_size = nstates2[3];
0173       
0174       unsigned int track2_silicon_cluster_size = std::numeric_limits<unsigned int>::quiet_NaN();
0175       if (siliconseed2)
0176       {
0177         track2_silicon_cluster_size = siliconseed2->size_cluster_keys();
0178       }
0179 
0180       std::vector<TrkrDefs::cluskey> ckeys2;
0181       if (siliconseed2)
0182       {
0183         ckeys2.insert(ckeys2.end(), siliconseed2->begin_cluster_keys(), siliconseed2->end_cluster_keys());
0184       }
0185 
0186       unsigned int track2_mvtx_cluster_size = 0;
0187       unsigned int track2_intt_cluster_size = 0;
0188       for (const auto& ckey : ckeys2)
0189       {
0190          auto detid = TrkrDefs::getTrkrId(ckey);
0191         if (detid == TrkrDefs::TrkrId::mvtxId)
0192         {
0193           track2_mvtx_cluster_size++;
0194         }
0195         else if (detid == TrkrDefs::TrkrId::inttId)
0196         {
0197           track2_intt_cluster_size++;
0198         }
0199       }
0200 
0201       // dca xy and dca z cut here compare to track dca cut
0202       Acts::Vector3 pos2(tr2->get_x(), tr2->get_y(), tr2->get_z());
0203       Acts::Vector3 mom2(tr2->get_px(), tr2->get_py(), tr2->get_pz());
0204       Acts::Vector3 dcaVals2 = calculateDca(tr2, mom2, pos2);
0205 
0206       if (fabs(dcaVals2(0)) < this_dca_cut2 or fabs(dcaVals2(1)) < this_dca_cut2)
0207       {
0208         continue;
0209       }
0210 
0211       // find DCA of these two tracks
0212       if (Verbosity() > 3)
0213       {
0214         std::cout << "Check DCA for tracks " << id1 << " and  " << id2 << std::endl;
0215       }
0216 
0217       if (tr1->get_charge() == tr2->get_charge())
0218       {
0219         // continue;
0220       }
0221 
0222       // declare these variables to pass into findPCAtwoTracks and fillHistogram by reference
0223       double pair_dca;
0224       Acts::Vector3 pca_rel1;
0225       Acts::Vector3 pca_rel2;
0226       double invariantMass;
0227       double invariantPt;
0228       float invariantPhi;
0229       float rapidity;
0230       float pseudorapidity;
0231 
0232       // Initial calculation of point of closest approach between the two tracks
0233       // This presently assumes straight line tracks to get a rough answer
0234       // Should update to use circles instead?
0235       findPcaTwoTracks(pos1, pos2, mom1, mom2, pca_rel1, pca_rel2, pair_dca);
0236 
0237       // tracks with small relative pca are k short candidates
0238       if (abs(pair_dca) < pair_dca_cut)
0239       {
0240         // Pair pca and dca were calculated with nominal track parameters and are approximate
0241         // Project tracks to this rough pca
0242         Eigen::Vector3d projected_pos1;
0243         Eigen::Vector3d projected_mom1;
0244         Eigen::Vector3d projected_pos2;
0245         Eigen::Vector3d projected_mom2;
0246 
0247         bool ret1 = projectTrackToPoint(tr1, pca_rel1, projected_pos1, projected_mom1);
0248         bool ret2 = projectTrackToPoint(tr2, pca_rel2, projected_pos2, projected_mom2);
0249 
0250         if (!ret1 or !ret2)
0251         {
0252           continue;
0253         }
0254 
0255         // recalculate pca starting with projected position and momentum
0256         double pair_dca_proj;
0257         Acts::Vector3 pca_rel1_proj;
0258         Acts::Vector3 pca_rel2_proj;
0259         findPcaTwoTracks(projected_pos1, projected_pos2, projected_mom1, projected_mom2, pca_rel1_proj, pca_rel2_proj, pair_dca_proj);
0260 
0261         // if(pair_dca_proj > pair_dca_cut) continue;
0262 
0263         // invariant mass is calculated in this method
0264         fillHistogram(projected_mom1, projected_mom2, recomass, invariantMass, invariantPt, invariantPhi, rapidity, pseudorapidity);
0265         fillNtp(tr1, tr2, dcaVals1, dcaVals2, pca_rel1, pca_rel2, pair_dca, invariantMass, invariantPt, invariantPhi, rapidity, pseudorapidity, projected_pos1, projected_pos2, projected_mom1, projected_mom2, pca_rel1_proj, pca_rel2_proj, pair_dca_proj, track1_silicon_cluster_size, track2_silicon_cluster_size, track1_mvtx_cluster_size, track1_mvtx_state_size, track1_intt_cluster_size, track1_intt_state_size, track2_mvtx_cluster_size, track2_mvtx_state_size, track2_intt_cluster_size, track2_intt_state_size, m_runNumber, m_evtNumber);
0266 
0267         if (Verbosity() > 1)
0268         {
0269           std::cout << " Accepted Track Pair" << std::endl;
0270       std::cout << " id1 " << id1 << " id2 " << id2 << std::endl;
0271       std::cout << " crossing1 " << crossing1 << " crossing2 " << crossing2 << std::endl;
0272           std::cout << " invariant mass: " << invariantMass << std::endl;
0273           std::cout << " track1 dca_cut: " << this_dca_cut << " track2 dca_cut: " << this_dca_cut2 << std::endl;
0274           std::cout << " dca3dxy1,dca3dz1,phi1: " << dcaVals1 << std::endl;
0275           std::cout << " dca3dxy2,dca3dz2,phi2: " << dcaVals2 << std::endl;
0276           std::cout << "Initial:  pca_rel1: " << pca_rel1 << " pca_rel2: " << pca_rel2 << std::endl;
0277           std::cout << " Initial: mom1: " << mom1 << " mom2: " << mom2 << std::endl;
0278           std::cout << "Proj_pca_rel:  proj_pos1: " << projected_pos1 << " proj_pos2: " << projected_pos2 << " proj_mom1: " << projected_mom1 << " proj_mom2: " << projected_mom2 << std::endl;
0279           std::cout << " Relative PCA = " << abs(pair_dca) << " pca_cut = " << pair_dca_cut << std::endl;
0280           std::cout << " charge 1: " << tr1->get_charge() << " charge2: " << tr2->get_charge() << std::endl;
0281           std::cout << "found viable projection" << std::endl;
0282           std::cout << "Final: pca_rel1_proj: " << pca_rel1_proj << " pca_rel2_proj: " << pca_rel2_proj << " mom1: " << projected_mom1 << " mom2: " << projected_mom2 << std::endl
0283                     << std::endl;
0284         }
0285 
0286         if (m_save_tracks)
0287         {
0288           m_output_trackMap = findNode::getClass<SvtxTrackMap>(topNode, m_output_trackMap_node_name);
0289           m_output_trackMap->insertWithKey(tr1, tr1->get_id());
0290           m_output_trackMap->insertWithKey(tr2, tr2->get_id());
0291         }
0292 
0293       }
0294     }
0295   }
0296   return 0;
0297 }
0298 
0299 std::vector<unsigned int> KshortReconstruction::getTrackStates(SvtxTrack *track)
0300 {
0301   std::vector<unsigned int> nstates;
0302   unsigned int nmapsstate = 0;
0303   unsigned int ninttstate = 0;
0304   unsigned int ntpcstate = 0;
0305   unsigned int nmmsstate = 0;
0306   
0307   // the track states from the Acts fit are fitted to fully corrected clusters, and are on the surface
0308   for (auto state_iter = track->begin_states();
0309        state_iter != track->end_states();
0310        ++state_iter)
0311     {
0312       SvtxTrackState* tstate = state_iter->second;
0313       auto stateckey = tstate->get_cluskey();
0314 
0315       switch (TrkrDefs::getTrkrId(stateckey))
0316     {
0317     case TrkrDefs::mvtxId:
0318       nmapsstate++;
0319       break;
0320     case TrkrDefs::inttId:
0321       ninttstate++;
0322       break;
0323     case TrkrDefs::tpcId:
0324       ntpcstate++;
0325       break;
0326     case TrkrDefs::micromegasId:
0327       nmmsstate++;
0328       break;
0329     default:
0330       std::cout << PHWHERE << " unknown key " << stateckey << std::endl;
0331       gSystem->Exit(1);
0332       exit(1);
0333     }
0334     }
0335   nstates.push_back(nmapsstate);
0336   nstates.push_back(ninttstate);
0337   nstates.push_back(ntpcstate);
0338   nstates.push_back(nmmsstate);
0339 
0340   return nstates;
0341 }
0342 
0343 void KshortReconstruction::fillNtp(SvtxTrack* track1, SvtxTrack* track2, Acts::Vector3 dcavals1, Acts::Vector3 dcavals2, Acts::Vector3 pca_rel1, Acts::Vector3 pca_rel2, double pair_dca, double invariantMass, double invariantPt, float invariantPhi, float rapidity, float pseudorapidity, Eigen::Vector3d projected_pos1, Eigen::Vector3d projected_pos2, Eigen::Vector3d projected_mom1, Eigen::Vector3d projected_mom2, Acts::Vector3 pca_rel1_proj, Acts::Vector3 pca_rel2_proj, double pair_dca_proj, unsigned int track1_silicon_cluster_size, unsigned int track2_silicon_cluster_size, unsigned int track1_mvtx_cluster_size,  unsigned int track1_mvtx_state_size, unsigned int track1_intt_cluster_size,  unsigned int track1_intt_state_size, unsigned int track2_mvtx_cluster_size,  unsigned int track2_mvtx_state_size, unsigned int track2_intt_cluster_size,  unsigned int track2_intt_state_size, int runNumber, int eventNumber)
0344 {
0345   double px1 = track1->get_px();
0346   double py1 = track1->get_py();
0347   double pz1 = track1->get_pz();
0348   auto *tpcSeed1 = track1->get_tpc_seed();
0349   size_t tpcClusters1 = tpcSeed1->size_cluster_keys();
0350   double eta1 = asinh(pz1 / sqrt(pow(px1, 2) + pow(py1, 2)));
0351 
0352   double px2 = track2->get_px();
0353   double py2 = track2->get_py();
0354   double pz2 = track2->get_pz();
0355   auto *tpcSeed2 = track2->get_tpc_seed();
0356   size_t tpcClusters2 = tpcSeed2->size_cluster_keys();
0357   double eta2 = asinh(pz2 / sqrt(pow(px2, 2) + pow(py2, 2)));
0358 
0359   auto vtxid = track1->get_vertex_id();
0360 
0361   Acts::Vector3 vertex(0, 0, track1->get_z());  // fake primary vertex
0362   auto *svtxVertex = m_vertexMap->get(vtxid);
0363   if (svtxVertex)
0364   {
0365     vertex(0) = svtxVertex->get_x();
0366     vertex(1) =  svtxVertex->get_y();
0367     vertex(2) = svtxVertex->get_z(); 
0368   }
0369 
0370   Acts::Vector3 pathLength = (pca_rel1 + pca_rel2) * 0.5 - vertex;
0371   Acts::Vector3 pathLength_proj = (pca_rel1_proj + pca_rel2_proj) * 0.5 - vertex;
0372 
0373   float mag_pathLength = sqrt(pow(pathLength(0), 2) + pow(pathLength(1), 2) + pow(pathLength(2), 2));
0374   float mag_pathLength_proj = sqrt(pow(pathLength_proj(0), 2) + pow(pathLength_proj(1), 2) + pow(pathLength_proj(2), 2));
0375 
0376   Acts::Vector3 projected_momentum = projected_mom1 + projected_mom2;
0377   float cos_theta_reco = pathLength_proj.dot(projected_momentum) / (projected_momentum.norm() * pathLength_proj.norm());
0378 
0379 
0380   float reco_info[] = {(float) track1->get_id(), (float) track1->get_crossing(), track1->get_x(), track1->get_y(), track1->get_z(), track1->get_px(), track1->get_py(), track1->get_pz(), (float) dcavals1(0), (float) dcavals1(1), (float) dcavals1(2), (float) pca_rel1(0), (float) pca_rel1(1), (float) pca_rel1(2), (float) eta1, (float) track1->get_charge(), (float) tpcClusters1, (float) track2->get_id(), (float) track2->get_crossing(), track2->get_x(), track2->get_y(), track2->get_z(), track2->get_px(), track2->get_py(), track2->get_pz(), (float) dcavals2(0), (float) dcavals2(1), (float) dcavals2(2), (float) pca_rel2(0), (float) pca_rel2(1), (float) pca_rel2(2), (float) eta2, (float) track2->get_charge(), (float) tpcClusters2, (float) vertex(0), (float) vertex(1), (float) vertex(2), (float) pair_dca, (float) invariantMass, (float) invariantPt, invariantPhi, (float) pathLength(0), (float) pathLength(1), (float) pathLength(2), mag_pathLength, rapidity, pseudorapidity, (float) projected_pos1(0), (float) projected_pos1(1), (float) projected_pos1(2), (float) projected_pos2(0), (float) projected_pos2(1), (float) projected_pos2(2), (float) projected_mom1(0), (float) projected_mom1(1), (float) projected_mom1(2), (float) projected_mom2(0), (float) projected_mom2(1), (float) projected_mom2(2), (float) pca_rel1_proj(0), (float) pca_rel1_proj(1), (float) pca_rel1_proj(2), (float) pca_rel2_proj(0), (float) pca_rel2_proj(1), (float) pca_rel2_proj(2), (float) pair_dca_proj, (float) pathLength_proj(0), (float) pathLength_proj(1), (float) pathLength_proj(2), mag_pathLength_proj, track1->get_quality(), track2->get_quality(), cos_theta_reco, (float) track1_silicon_cluster_size, (float) track2_silicon_cluster_size, (float) track1_mvtx_cluster_size, (float) track1_mvtx_state_size, (float) track1_intt_cluster_size,  (float) track1_intt_state_size, (float) track2_mvtx_cluster_size, (float) track2_mvtx_state_size,  (float) track2_intt_cluster_size, (float) track2_intt_state_size, (float) runNumber, (float) eventNumber};
0381 
0382   ntp_reco_info->Fill(reco_info);
0383 }
0384 
0385 void KshortReconstruction::fillHistogram(Eigen::Vector3d mom1, Eigen::Vector3d mom2, TH1* massreco, double& invariantMass, double& invariantPt, float& invariantPhi, float& rapidity, float& pseudorapidity)
0386 {
0387   double E1 = sqrt(pow(mom1(0), 2) + pow(mom1(1), 2) + pow(mom1(2), 2) + pow(decaymass, 2));
0388   double E2 = sqrt(pow(mom2(0), 2) + pow(mom2(1), 2) + pow(mom2(2), 2) + pow(decaymass, 2));
0389 
0390   TLorentzVector v1(mom1(0), mom1(1), mom1(2), E1);
0391   TLorentzVector v2(mom2(0), mom2(1), mom2(2), E2);
0392 
0393   TLorentzVector tsum;
0394   tsum = v1 + v2;
0395 
0396   rapidity = tsum.Rapidity();
0397   pseudorapidity = tsum.Eta();
0398   invariantMass = tsum.M();
0399   invariantPt = tsum.Pt();
0400   invariantPhi = tsum.Phi();
0401 
0402   if (Verbosity() > 2)
0403   {
0404     std::cout << "px1: " << mom1(0) << " py1: " << mom1(1) << " pz1: " << mom1(2) << " E1: " << E1 << std::endl;
0405     std::cout << "px2: " << mom2(0) << " py2: " << mom2(1) << " pz2: " << mom2(2) << " E2: " << E2 << std::endl;
0406     std::cout << "tsum: " << tsum(0) << " " << tsum(1) << " " << tsum(2) << " " << tsum(3) << std::endl;
0407     std::cout << "invariant mass: " << invariantMass << " invariant Pt: " << invariantPt << " invariantPhi: " << invariantPhi << std::endl;
0408   }
0409 
0410   if (invariantPt > invariant_pt_cut)
0411   {
0412     massreco->Fill(invariantMass);
0413   }
0414 }
0415 
0416 bool KshortReconstruction::projectTrackToPoint(SvtxTrack* track, Eigen::Vector3d PCA, Eigen::Vector3d& pos, Eigen::Vector3d& mom)
0417 {
0418   bool ret = true;
0419 
0420   /// create perigee surface
0421   ActsPropagator actsPropagator(_tGeometry);
0422   auto perigee = actsPropagator.makeVertexSurface(PCA);  // PCA is in cm here
0423   auto params = actsPropagator.makeTrackParams(track, m_vertexMap);
0424   if (!params.ok())
0425   {
0426     return false;
0427   }
0428   auto result = actsPropagator.propagateTrack(params.value(), perigee);
0429 
0430   if (result.ok())
0431   {
0432     auto projectionPos = result.value().second.position(_tGeometry->geometry().getGeoContext());
0433     const auto momentum = result.value().second.momentum();
0434     pos(0) = projectionPos.x() / Acts::UnitConstants::cm;
0435     pos(1) = projectionPos.y() / Acts::UnitConstants::cm;
0436     pos(2) = projectionPos.z() / Acts::UnitConstants::cm;
0437 
0438     if (Verbosity() > 2)
0439     {
0440       std::cout << "                 Input PCA " << PCA << "  projection out " << pos << std::endl;
0441     }
0442 
0443     mom(0) = momentum.x();
0444     mom(1) = momentum.y();
0445     mom(2) = momentum.z();
0446   }
0447   else
0448   {
0449     pos(0) = track->get_x();
0450     pos(1) = track->get_y();
0451     pos(2) = track->get_z();
0452 
0453     mom(0) = track->get_px();
0454     mom(1) = track->get_py();
0455     mom(2) = track->get_pz();
0456 
0457     if(Verbosity() > 0)
0458       {
0459     std::cout << result.error() << std::endl;
0460     std::cout << result.error().message() << std::endl;
0461     std::cout << " Failed projection of track with: " << std::endl;
0462     std::cout << " x,y,z = " << track->get_x() << "  " << track->get_y() << "  " << track->get_z() << std::endl;
0463     std::cout << " px,py,pz = " << track->get_px() << "  " << track->get_py() << "  " << track->get_pz() << std::endl;
0464     std::cout << " to point (x,y,z) = " << PCA(0) / Acts::UnitConstants::cm << "  " << PCA(1) / Acts::UnitConstants::cm << "  " << PCA(2) / Acts::UnitConstants::cm << std::endl;
0465       }
0466 
0467     //    ret = false;
0468   }
0469 
0470   return ret;
0471 }
0472 
0473 bool KshortReconstruction::projectTrackToCylinder(SvtxTrack* track, double Radius, Eigen::Vector3d& pos, Eigen::Vector3d& mom)
0474 {
0475   // Make a cylinder surface at the radius and project the track to that
0476   bool ret = true;
0477   const double eta = 2.0;
0478   const double theta = 2. * atan(exp(-eta));
0479   const double halfZ = Radius / tan(theta) * Acts::UnitConstants::cm;
0480   Radius *= Acts::UnitConstants::cm;
0481 
0482   /// Make a cylindrical surface at (0,0,0) aligned along the z axis
0483   auto transform = Acts::Transform3::Identity();
0484 
0485   std::shared_ptr<Acts::CylinderSurface> cylSurf =
0486       Acts::Surface::makeShared<Acts::CylinderSurface>(transform,
0487                                                        Radius,
0488                                                        halfZ);
0489   ActsPropagator actsPropagator(_tGeometry);
0490   auto params = actsPropagator.makeTrackParams(track, m_vertexMap);
0491   if (!params.ok())
0492   {
0493     return false;
0494   }
0495 
0496   auto result = actsPropagator.propagateTrack(params.value(), cylSurf);
0497   if (result.ok())
0498   {
0499     auto projectionPos = result.value().second.position(_tGeometry->geometry().getGeoContext());
0500     const auto momentum = result.value().second.momentum();
0501     pos(0) = projectionPos.x() / Acts::UnitConstants::cm;
0502     pos(1) = projectionPos.y() / Acts::UnitConstants::cm;
0503     pos(2) = projectionPos.z() / Acts::UnitConstants::cm;
0504 
0505     mom(0) = momentum.x();
0506     mom(1) = momentum.y();
0507     mom(2) = momentum.z();
0508   }
0509   else
0510   {
0511     ret = false;
0512   }
0513 
0514   return ret;
0515 }
0516 
0517 Acts::Vector3 KshortReconstruction::getVertex(SvtxTrack* track)
0518 {
0519   auto vertexId = track->get_vertex_id();
0520   const SvtxVertex* svtxVertex = m_vertexMap->get(vertexId);
0521   Acts::Vector3 vertex = Acts::Vector3::Zero();
0522   if (svtxVertex)
0523   {
0524     vertex(0) = svtxVertex->get_x() * Acts::UnitConstants::cm;
0525     vertex(1) = svtxVertex->get_y() * Acts::UnitConstants::cm;
0526     vertex(2) = svtxVertex->get_z() * Acts::UnitConstants::cm;
0527   }
0528 
0529   return vertex;
0530 }
0531 
0532 void KshortReconstruction::findPcaTwoTracks(const Acts::Vector3& pos1, const Acts::Vector3& pos2, Acts::Vector3 mom1, Acts::Vector3 mom2, Acts::Vector3& pca1, Acts::Vector3& pca2, double& dca) const
0533 {
0534   TLorentzVector v1;
0535   TLorentzVector v2;
0536 
0537   double px1 = mom1(0);
0538   double py1 = mom1(1);
0539   double pz1 = mom1(2);
0540   double px2 = mom2(0);
0541   double py2 = mom2(1);
0542   double pz2 = mom2(2);
0543 
0544   Float_t E1 = sqrt(pow(px1, 2) + pow(py1, 2) + pow(pz1, 2) + pow(decaymass, 2));
0545   Float_t E2 = sqrt(pow(px2, 2) + pow(py2, 2) + pow(pz2, 2) + pow(decaymass, 2));
0546 
0547   v1.SetPxPyPzE(px1, py1, pz1, E1);
0548   v2.SetPxPyPzE(px2, py2, pz2, E2);
0549 
0550   // calculate lorentz vector
0551   const Eigen::Vector3d& a1 = pos1;
0552   const Eigen::Vector3d& a2 = pos2;
0553 
0554   Eigen::Vector3d b1(v1.Px(), v1.Py(), v1.Pz());
0555   Eigen::Vector3d b2(v2.Px(), v2.Py(), v2.Pz());
0556 
0557   // The shortest distance between two skew lines described by
0558   //  a1 + c * b1
0559   //  a2 + d * b2
0560   // where a1, a2, are vectors representing points on the lines, b1, b2 are direction vectors, and c and d are scalars
0561   // dca = (b1 x b2) .(a2-a1) / |b1 x b2|
0562 
0563   // bcrossb/mag_bcrossb is a unit vector perpendicular to both direction vectors b1 and b2
0564   auto bcrossb = b1.cross(b2);
0565   auto mag_bcrossb = bcrossb.norm();
0566   // a2-a1 is the vector joining any arbitrary points on the two lines
0567   auto aminusa = a2 - a1;
0568 
0569   // The DCA of these two lines is the projection of a2-a1 along the direction of the perpendicular to both
0570   // remember that a2-a1 is longer than (or equal to) the dca by definition
0571   dca = 999;
0572   if (mag_bcrossb != 0)
0573   {
0574     dca = bcrossb.dot(aminusa) / mag_bcrossb;
0575   }
0576   else
0577   {
0578     return;  // same track, skip combination
0579   }
0580 
0581   // get the points at which the normal to the lines intersect the lines, where the lines are perpendicular
0582   double X = b1.dot(b2) - (b1.dot(b1) * b2.dot(b2) / b2.dot(b1));
0583   double Y = (a2.dot(b2) - a1.dot(b2)) - ((a2.dot(b1) - a1.dot(b1)) * b2.dot(b2) / b2.dot(b1));
0584   double c = Y / X;
0585 
0586   double F = b1.dot(b1) / b2.dot(b1);
0587   double G = -(a2.dot(b1) - a1.dot(b1)) / b2.dot(b1);
0588   double d = (c * F) + G;
0589 
0590   // then the points of closest approach are:
0591   pca1 = a1 + c * b1;
0592   pca2 = a2 + d * b2;
0593 
0594   return;
0595 }
0596 
0597 KshortReconstruction::KshortReconstruction(const std::string& name)
0598   : SubsysReco(name)
0599 {
0600 }
0601 
0602 Acts::Vector3 KshortReconstruction::calculateDca(SvtxTrack* track, const Acts::Vector3& momentum, Acts::Vector3 position)
0603 {
0604   // For the purposes of this module, we set default values to prevent this track from being rejected if the dca calc fails
0605   Acts::Vector3 r = momentum.cross(Acts::Vector3(0., 0., 1.));
0606   float phi = atan2(r(1), r(0));
0607   Acts::Vector3 outVals(track_dca_cut*1.1, track_dca_cut*1.1, phi);
0608   auto vtxid = track->get_vertex_id();
0609   if (!m_vertexMap)
0610   {
0611     //std::cout << "Could not find m_vertexmap " << std::endl;
0612     return outVals;
0613   }
0614   auto *svtxVertex = m_vertexMap->get(vtxid);
0615   if (!svtxVertex)
0616   {
0617     //std::cout << "Could not find vtxid in m_vertexMap " << vtxid << std::endl;
0618     return outVals;
0619   }
0620   Acts::Vector3 vertex(svtxVertex->get_x(), svtxVertex->get_y(), svtxVertex->get_z());
0621   position -= vertex;
0622 
0623   Acts::RotationMatrix3 rot;
0624   rot(0, 0) = std::cos(phi);
0625   rot(0, 1) = -std::sin(phi);
0626   rot(0, 2) = 0;
0627   rot(1, 0) = std::sin(phi);
0628   rot(1, 1) = std::cos(phi);
0629   rot(1, 2) = 0;
0630   rot(2, 0) = 0;
0631   rot(2, 1) = 0;
0632   rot(2, 2) = 1;
0633 
0634   Acts::Vector3 pos_R = rot * position;
0635   double dca3dxy = pos_R(0);
0636   double dca3dz = pos_R(2);
0637 
0638   outVals(0) = abs(dca3dxy);
0639   outVals(1) = abs(dca3dz);
0640   outVals(2) = phi;
0641 
0642   if (Verbosity() > 4)
0643   {
0644     std::cout << " pre-position: " << position << std::endl;
0645     std::cout << " vertex: " << vertex << std::endl;
0646     std::cout << " vertex subtracted-position: " << position << std::endl;
0647   }
0648 
0649   return outVals;
0650 }
0651 
0652 int KshortReconstruction::InitRun(PHCompositeNode* topNode)
0653 {
0654   const char* cfilepath = filepath.c_str();
0655   fout = new TFile(cfilepath, "recreate");
0656   ntp_reco_info = new TNtuple("ntp_reco_info", "decay_pairs", "id1:crossing1:x1:y1:z1:px1:py1:pz1:dca3dxy1:dca3dz1:phi1:pca_rel1_x:pca_rel1_y:pca_rel1_z:eta1:charge1:tpcClusters_1:id2:crossing2:x2:y2:z2:px2:py2:pz2:dca3dxy2:dca3dz2:phi2:pca_rel2_x:pca_rel2_y:pca_rel2_z:eta2:charge2:tpcClusters_2:vertex_x:vertex_y:vertex_z:pair_dca:invariant_mass:invariant_pt:invariantPhi:pathlength_x:pathlength_y:pathlength_z:pathlength:rapidity:pseudorapidity:projected_pos1_x:projected_pos1_y:projected_pos1_z:projected_pos2_x:projected_pos2_y:projected_pos2_z:projected_mom1_x:projected_mom1_y:projected_mom1_z:projected_mom2_x:projected_mom2_y:projected_mom2_z:projected_pca_rel1_x:projected_pca_rel1_y:projected_pca_rel1_z:projected_pca_rel2_x:projected_pca_rel2_y:projected_pca_rel2_z:projected_pair_dca:projected_pathlength_x:projected_pathlength_y:projected_pathlength_z:projected_pathlength:quality1:quality2:cosThetaReco:track1_silicon_clusters:track2_silicon_clusters:track1_mvtx_clusters:track1_mvtx_states:track1_intt_clusters:track1_intt_states:track2_mvtx_clusters:track2_mvtx_states:track2_intt_clusters:track2_intt_states:runNumber:eventNumber");
0657 
0658   getNodes(topNode);
0659 
0660   recomass = new TH1D("recomass", "recomass", 1000, 0.0, 1);  // root histogram arguments: name,title,bins,minvalx,maxvalx
0661 
0662   //Add new track map to save selected tracks
0663   if (m_save_tracks)
0664   {
0665     PHNodeIterator nodeIter(topNode);
0666 
0667     PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(nodeIter.findFirst("PHCompositeNode", "DST"));
0668     if (!dstNode)
0669     {
0670       dstNode = new PHCompositeNode("DST");
0671       topNode->addNode(dstNode);
0672       std::cout << "DST node added" << std::endl;
0673     }
0674 
0675     m_output_trackMap = new SvtxTrackMap_v2();
0676     PHIODataNode<PHObject> *outputTrackNode = new PHIODataNode<PHObject>(m_output_trackMap, m_output_trackMap_node_name, "PHObject");
0677     dstNode->addNode(outputTrackNode);
0678     if (Verbosity() > 1) { std::cout << m_output_trackMap_node_name << " node added" << std::endl; }
0679   }
0680 
0681   return 0;
0682 }
0683 
0684 int KshortReconstruction::End(PHCompositeNode* /**topNode*/)
0685 {
0686   fout->cd();
0687   ntp_reco_info->Write();
0688   recomass->Write();
0689   fout->Close();
0690 
0691   return 0;
0692 }
0693 
0694 int KshortReconstruction::getNodes(PHCompositeNode* topNode)
0695 {
0696   m_svtxTrackMap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0697   if (!m_svtxTrackMap)
0698   {
0699     std::cout << PHWHERE << "No SvtxTrackMap on node tree, exiting." << std::endl;
0700     return Fun4AllReturnCodes::ABORTEVENT;
0701   }
0702 
0703   m_vertexMap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0704   if (!m_vertexMap)
0705   {
0706     std::cout << PHWHERE << "No vertexMap on node tree, exiting." << std::endl;
0707     return Fun4AllReturnCodes::ABORTEVENT;
0708   }
0709 
0710   _tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0711   if (!_tGeometry)
0712   {
0713     std::cout << PHWHERE << "Error, can't find acts tracking geometry" << std::endl;
0714     return Fun4AllReturnCodes::ABORTEVENT;
0715   }
0716 
0717   return Fun4AllReturnCodes::EVENT_OK;
0718 }
[ 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!