packages/trackreco/PHCosmicsTrkFitter.cc

0001 #include "PHCosmicsTrkFitter.h"
0002 #include "MakeSourceLinks.h"
0003
0004 /// Tracking includes
0005 #include <trackbase/Calibrator.h>
0006 #include <trackbase/ClusterErrorPara.h>
0007 #include <trackbase/InttDefs.h>
0008 #include <trackbase/MvtxDefs.h>
0009 #include <trackbase/TpcDefs.h>
0010 #include <trackbase/TrackFitUtils.h>
0011 #include <trackbase/TrkrCluster.h>
0012 #include <trackbase/TrkrClusterContainer.h>
0013
0014 #include <trackbase_historic/ActsTransformations.h>
0015 #include <trackbase_historic/SvtxAlignmentStateMap_v1.h>
0016 #include <trackbase_historic/SvtxTrackMap_v2.h>
0017 #include <trackbase_historic/SvtxTrackState_v1.h>
0018 #include <trackbase_historic/SvtxTrack_v4.h>
0019 #include <trackbase_historic/TrackSeed.h>
0020 #include <trackbase_historic/TrackSeedContainer.h>
0021 #include <trackbase_historic/TrackSeedHelper.h>
0022
0023 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0024
0025 #include <micromegas/MicromegasDefs.h>
0026
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028 #include <phool/PHCompositeNode.h>
0029 #include <phool/PHDataNode.h>
0030 #include <phool/PHNode.h>
0031 #include <phool/PHNodeIterator.h>
0032 #include <phool/PHObject.h>
0033 #include <phool/PHTimer.h>
0034 #include <phool/getClass.h>
0035 #include <phool/phool.h>
0036
0037 #include <tpc/TpcDistortionCorrectionContainer.h>
0038
0039 #include <Acts/EventData/MultiTrajectory.hpp>
0040 #include <Acts/EventData/MultiTrajectoryHelpers.hpp>
0041 #include <Acts/EventData/TrackParameters.hpp>
0042 #include <Acts/Surfaces/PerigeeSurface.hpp>
0043 #include <Acts/Surfaces/PlaneSurface.hpp>
0044 #include <Acts/Surfaces/Surface.hpp>
0045 #include <Acts/TrackFitting/GainMatrixSmoother.hpp>
0046 #include <Acts/TrackFitting/GainMatrixUpdater.hpp>
0047
0048 #include <TTree.h>
0049 #include <TFile.h>
0050 #include <TVector3.h>
0051
0052 #include <cmath>
0053 #include <iostream>
0054 #include <vector>
0055
0056 namespace
0057 {
0058   // check vector validity
0059   inline bool is_valid(const Acts::Vector3& vec)
0060   {
0061     return !(std::isnan(vec.x()) || std::isnan(vec.y()) || std::isnan(vec.z()));
0062   }
0063   template <class T>
0064   inline T square(const T& x)
0065   {
0066     return x * x;
0067   }
0068
0069   /// get radius from coordinates
0070   template <class T>
0071   T radius(const T& x, const T& y)
0072   {
0073     return std::sqrt(square(x) + square(y));
0074   }
0075
0076 }  // namespace
0077
0078 #include <trackbase/alignmentTransformationContainer.h>
0079 #include <Eigen/Dense>
0080 #include <Eigen/Geometry>
0081
0082 PHCosmicsTrkFitter::PHCosmicsTrkFitter(const std::string& name)
0083   : SubsysReco(name)
0084   , m_trajectories(nullptr)
0085 {
0086 }
0087
0088 int PHCosmicsTrkFitter::InitRun(PHCompositeNode* topNode)
0089 {
0090   if (Verbosity() > 1)
0091   {
0092     std::cout << "Setup PHCosmicsTrkFitter" << std::endl;
0093   }
0094
0095   if (createNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0096   {
0097     return Fun4AllReturnCodes::ABORTEVENT;
0098   }
0099
0100   if (getNodes(topNode) != Fun4AllReturnCodes::EVENT_OK)
0101   {
0102     return Fun4AllReturnCodes::ABORTEVENT;
0103   }
0104
0105   // configure alignStates
0106   m_alignStates.loadNodes(topNode);
0107   m_alignStates.verbosity(Verbosity());
0108   m_alignStates.fieldMap(m_fieldMap);
0109
0110   // detect const field
0111   std::istringstream stringline(m_fieldMap);
0112   stringline >> fieldstrength;
0113   if (!stringline.fail())  // it is a float
0114   {
0115     m_ConstField = true;
0116   }
0117
0118   m_fitCfg.fit = ActsTrackFittingAlgorithm::makeKalmanFitterFunction(
0119       m_tGeometry->geometry().tGeometry,
0120       m_tGeometry->geometry().magField);
0121
0122   m_outlierFinder.verbosity = Verbosity();
0123   std::map<long unsigned int, float> chi2Cuts;
0124   chi2Cuts.insert(std::make_pair(10, 4));
0125   chi2Cuts.insert(std::make_pair(12, 4));
0126   chi2Cuts.insert(std::make_pair(14, 9));
0127   chi2Cuts.insert(std::make_pair(16, 4));
0128   m_outlierFinder.chi2Cuts = chi2Cuts;
0129   if (m_useOutlierFinder)
0130   {
0131     m_fitCfg.fit->outlierFinder(m_outlierFinder);
0132   }
0133
0134   _tpccellgeo = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0135
0136   if (m_actsEvaluator)
0137   {
0138     m_evaluator = std::make_unique<ActsEvaluator>(m_evalname);
0139     m_evaluator->Init(topNode);
0140     m_evaluator->verbosity(Verbosity());
0141   }
0142
0143   if (m_seedClusAnalysis)
0144   {
0145     m_outfile = new TFile(m_evalname.c_str(), "RECREATE");
0146     m_tree = new TTree("seedclustree", "Tree with cosmic seeds and their clusters");
0147     makeBranches();
0148   }
0149
0150   if (Verbosity() > 1)
0151   {
0152     std::cout << "Finish PHCosmicsTrkFitter Setup" << std::endl;
0153   }
0154
0155   return Fun4AllReturnCodes::EVENT_OK;
0156 }
0157
0158 int PHCosmicsTrkFitter::process_event(PHCompositeNode* topNode)
0159 {
0160   auto logLevel = Acts::Logging::FATAL;
0161
0162   if (m_actsEvaluator)
0163   {
0164     m_evaluator->next_event(topNode);
0165   }
0166
0167   if (Verbosity() > 1)
0168   {
0169     std::cout << PHWHERE << "Events processed: " << m_event << std::endl;
0170     std::cout << "Start PHCosmicsTrkFitter::process_event" << std::endl;
0171     if (Verbosity() > 30)
0172     {
0173       logLevel = Acts::Logging::VERBOSE;
0174     }
0175   }
0176
0177   /// Fill an additional track map if using the acts evaluator
0178   /// for proto track comparison to fitted track
0179   if (m_actsEvaluator)
0180   {
0181     /// wipe at the beginning of every new fit pass, so that the seeds
0182     /// are whatever is currently in SvtxTrackMap
0183     m_seedTracks->clear();
0184     for (const auto& [key, track] : *m_trackMap)
0185     {
0186       m_seedTracks->insert(track);
0187     }
0188   }
0189
0190   loopTracks(logLevel);
0191
0192   // put this in the output file
0193   if (Verbosity() > 0)
0194   {
0195     std::cout << " SvtxTrackMap size is now " << m_trackMap->size()
0196               << std::endl;
0197   }
0198
0199   m_event++;
0200   return Fun4AllReturnCodes::EVENT_OK;
0201 }
0202
0203 int PHCosmicsTrkFitter::ResetEvent(PHCompositeNode* /*topNode*/)
0204 {
0205   if (Verbosity() > 1)
0206   {
0207     std::cout << "Reset PHCosmicsTrkFitter" << std::endl;
0208   }
0209
0210   m_trajectories->clear();
0211
0212   return Fun4AllReturnCodes::EVENT_OK;
0213 }
0214
0215 int PHCosmicsTrkFitter::End(PHCompositeNode* /*topNode*/)
0216 {
0217   if (m_actsEvaluator)
0218   {
0219     m_evaluator->End();
0220   }
0221   if (m_seedClusAnalysis)
0222   {
0223     m_outfile->cd();
0224     m_tree->Write();
0225     m_outfile->Close();
0226   }
0227   if (Verbosity() > 0)
0228   {
0229     std::cout << "The Acts track fitter succeeded " << m_nGoodFits
0230               << " times and had " << m_nBadFits
0231               << " fits return an error" << std::endl;
0232     std::cout << "There were " << m_nLongSeeds << " long seeds" << std::endl;
0233     std::cout << "Finished PHCosmicsTrkFitter" << std::endl;
0234   }
0235   return Fun4AllReturnCodes::EVENT_OK;
0236 }
0237
0238 void PHCosmicsTrkFitter::loopTracks(Acts::Logging::Level logLevel)
0239 {
0240   auto logger = Acts::getDefaultLogger("PHCosmicsTrkFitter", logLevel);
0241
0242   if (Verbosity() > 0)
0243   {
0244     std::cout << " seed map size " << m_seedMap->size() << std::endl;
0245   }
0246
0247   for (auto track : *m_seedMap)
0248   {
0249     if (!track)
0250     {
0251       continue;
0252     }
0253
0254     unsigned int tpcid = track->get_tpc_seed_index();
0255     unsigned int siid = track->get_silicon_seed_index();
0256
0257     // get the crossing number
0258     auto siseed = m_siliconSeeds->get(siid);
0259     short crossing = 0;
0260
0261     auto tpcseed = m_tpcSeeds->get(tpcid);
0262     if (Verbosity() > 1)
0263     {
0264       std::cout << "TPC id " << tpcid << std::endl;
0265       std::cout << "Silicon id " << siid << std::endl;
0266     }
0267
0268     /// Need to also check that the tpc seed wasn't removed by the ghost finder
0269     if (!tpcseed)
0270     {
0271       continue;
0272     }
0273
0274     if (Verbosity() > 1)
0275     {
0276       const auto position = TrackSeedHelper::get_xyz(tpcseed);
0277       std::cout << " tpc seed position is (x,y,z) = " << position.x() << "  " << position.y() << "  " << position.z() << std::endl;
0278     }
0279
0280     ActsTrackFittingAlgorithm::MeasurementContainer measurements;
0281
0282     SourceLinkVec sourceLinks;
0283
0284     MakeSourceLinks makeSourceLinks;
0285     makeSourceLinks.initialize(_tpccellgeo);
0286     makeSourceLinks.setVerbosity(Verbosity());
0287     makeSourceLinks.set_pp_mode(false);
0288
0289     makeSourceLinks.resetTransientTransformMap(
0290         m_alignmentTransformationMapTransient,
0291         m_transient_id_set,
0292         m_tGeometry);
0293
0294     if (siseed)
0295     {
0296       // silicon source links
0297       sourceLinks = makeSourceLinks.getSourceLinks(
0298         siseed,
0299         measurements,
0300         m_clusterContainer,
0301         m_tGeometry,
0302         m_globalPositionWrapper,
0303         m_alignmentTransformationMapTransient,
0304         m_transient_id_set,
0305         crossing);
0306     }
0307
0308     // tpc source links
0309     const auto tpcSourceLinks = makeSourceLinks.getSourceLinks(
0310       tpcseed,
0311       measurements,
0312       m_clusterContainer,
0313       m_tGeometry,
0314       m_globalPositionWrapper,
0315       m_alignmentTransformationMapTransient,
0316       m_transient_id_set,
0317       crossing);
0318
0319     sourceLinks.insert(sourceLinks.end(), tpcSourceLinks.begin(), tpcSourceLinks.end());
0320
0321     if (sourceLinks.size() < 5)
0322     {
0323       continue;
0324     }
0325     int charge = 0;
0326     float cosmicslope = 0;
0327
0328     getCharge(tpcseed, charge, cosmicslope);
0329
0330     // copy transient map for this track into transient geoContext
0331     m_transient_geocontext = m_alignmentTransformationMapTransient;
0332
0333     {
0334       // get positions from cluster keys
0335       // TODO: should implement distortions
0336       TrackSeedHelper::position_map_t positions;
0337       for( auto key_iter = tpcseed->begin_cluster_keys(); key_iter != tpcseed->end_cluster_keys(); ++key_iter )
0338       {
0339         const auto& key(*key_iter);
0340         positions.emplace(key, m_tGeometry->getGlobalPosition( key, m_clusterContainer->findCluster(key)));
0341       }
0342
0343       TrackSeedHelper::circleFitByTaubin(tpcseed, positions, 0, 58);
0344     }
0345
0346     float tpcR = fabs(1. / tpcseed->get_qOverR());
0347     float tpcx = tpcseed->get_X0();
0348     float tpcy = tpcseed->get_Y0();
0349
0350     const auto intersect =
0351         TrackFitUtils::circle_circle_intersection(m_vertexRadius,
0352                                                   tpcR, tpcx, tpcy);
0353     float intx, inty;
0354
0355     if (std::get<1>(intersect) > std::get<3>(intersect))
0356     {
0357       intx = std::get<0>(intersect);
0358       inty = std::get<1>(intersect);
0359     }
0360     else
0361     {
0362       intx = std::get<2>(intersect);
0363       inty = std::get<3>(intersect);
0364     }
0365     std::vector<TrkrDefs::cluskey> keys;
0366     std::vector<Acts::Vector3> clusPos;
0367     std::copy(tpcseed->begin_cluster_keys(), tpcseed->end_cluster_keys(), std::back_inserter(keys));
0368     TrackFitUtils::getTrackletClusters(m_tGeometry, m_clusterContainer,
0369                                        clusPos, keys);
0370     TrackFitUtils::position_vector_t xypoints, rzpoints;
0371     for (auto& pos : clusPos)
0372     {
0373       float clusr = radius(pos.x(), pos.y());
0374       if (pos.y() < 0)
0375       {
0376         clusr *= -1;
0377       }
0378
0379       // exclude silicon and tpot clusters for now
0380       if (fabs(clusr) > 80 || fabs(clusr) < 30)
0381       {
0382         continue;
0383       }
0384       xypoints.push_back(std::make_pair(pos.x(), pos.y()));
0385       rzpoints.push_back(std::make_pair(pos.z(), clusr));
0386     }
0387
0388     auto rzparams = TrackFitUtils::line_fit(rzpoints);
0389     float fulllineintz = std::get<1>(rzparams);
0390     float fulllineslope = std::get<0>(rzparams);
0391
0392     float slope = tpcseed->get_slope();
0393     float intz = m_vertexRadius * slope + tpcseed->get_Z0();
0394
0395     Acts::Vector3 inter(intx, inty, intz);
0396
0397     std::vector<float> tpcparams{tpcR, tpcx, tpcy, tpcseed->get_slope(),
0398                                  tpcseed->get_Z0()};
0399     auto tangent = TrackFitUtils::get_helix_tangent(tpcparams,
0400                                                     inter);
0401
0402     auto tan = tangent.second;
0403     auto pca = tangent.first;
0404
0405     float p;
0406     if (m_ConstField)
0407     {
0408       p = cosh(tpcseed->get_eta()) * fabs(1. / tpcseed->get_qOverR()) * (0.3 / 100) * fieldstrength;
0409     }
0410     else
0411     {
0412       p = tpcseed->get_p();
0413     }
0414
0415     tan *= p;
0416
0417     //! if we got the opposite seed then z will be backwards, so we take the
0418     //! value of tan.z() multiplied by the sign of the slope determined for
0419     //! the full cosmic track
0420     //! same with px/py since a single cosmic produces two seeds that bend
0421     //! in opposite directions
0422     float theta = std::atan(fulllineslope);
0423     /// Normalize to 0<theta<pi
0424     if (theta < 0)
0425     {
0426       theta += M_PI;
0427     }
0428     float pz = std::cos(theta) * p;
0429     if (fulllineslope < 0)
0430     {
0431       pz = fabs(pz);
0432     }
0433     else
0434     {
0435       pz = fabs(pz) * -1;
0436     }
0437     Acts::Vector3 momentum = Acts::Vector3::Zero();
0438
0439     if (!m_zeroField)
0440     {
0441       momentum.x() = charge < 0 ? tan.x() : tan.x() * -1;
0442       momentum.y() = charge < 0 ? tan.y() : tan.y() * -1;
0443     }
0444     else
0445     {
0446       auto xyparams = TrackFitUtils::line_fit(xypoints);
0447       float fulllineslopexy = std::get<0>(xyparams);
0448       if (fulllineslopexy < 0)
0449       {
0450         momentum.x() = fabs(tan.x());
0451       }
0452       else
0453       {
0454         momentum.x() = fabs(tan.x()) * -1;
0455       }
0456       momentum.y() = fabs(tan.y()) * -1;
0457     }
0458
0459     momentum.z() = pz;
0460     Acts::Vector3 position(pca.x(), pca.y(),
0461                            (m_vertexRadius - fulllineintz) / fulllineslope);
0462
0463     position *= Acts::UnitConstants::cm;
0464     if (!is_valid(momentum))
0465     {
0466       continue;
0467     }
0468
0469     auto pSurface = Acts::Surface::makeShared<Acts::PerigeeSurface>(
0470         position);
0471     auto actsFourPos = Acts::Vector4(position(0), position(1),
0472                                      position(2),
0473                                      10 * Acts::UnitConstants::ns);
0474
0475     if (sourceLinks.size() > 20)
0476     {
0477       m_nLongSeeds++;
0478     }
0479     Acts::BoundSquareMatrix cov = setDefaultCovariance();
0480     if (m_seedClusAnalysis)
0481     {
0482       clearVectors();
0483       m_seed = tpcid;
0484       m_R = tpcR;
0485       m_X0 = tpcx;
0486       m_Y0 = tpcy;
0487       m_Z0 = fulllineintz;
0488       m_slope = fulllineslope;
0489       m_pcax = position(0);
0490       m_pcay = position(1);
0491       m_pcaz = position(2);
0492       m_px = momentum(0);
0493       m_py = momentum(1);
0494       m_pz = momentum(2);
0495       m_charge = charge;
0496       fillVectors(siseed, tpcseed);
0497       m_tree->Fill();
0498     }
0499     //! Reset the track seed with the dummy covariance
0500     auto seed = ActsTrackFittingAlgorithm::TrackParameters::create(
0501         pSurface,
0502         m_transient_geocontext,
0503         actsFourPos,
0504         momentum,
0505         charge / momentum.norm(),
0506         cov,
0507         Acts::ParticleHypothesis::muon(),
0508         100 * Acts::UnitConstants::cm);
0509     if (!seed.ok())
0510     {
0511       std::cout << "Could not create track params, skipping track" << std::endl;
0512       continue;
0513     }
0514
0515     if (Verbosity() > 0)
0516     {
0517       std::cout << "Source link size " << sourceLinks.size() << std::endl;
0518       printTrackSeed(seed.value());
0519     }
0520
0521     //! Set host of propagator options for Acts to do e.g. material integration
0522     Acts::PropagatorPlainOptions ppPlainOptions;
0523
0524     auto calibptr = std::make_unique<Calibrator>();
0525     CalibratorAdapter calibrator{*calibptr, measurements};
0526
0527     auto magcontext = m_tGeometry->geometry().magFieldContext;
0528     auto calibcontext = m_tGeometry->geometry().calibContext;
0529
0530     ActsTrackFittingAlgorithm::GeneralFitterOptions
0531         kfOptions{
0532             m_transient_geocontext,
0533             magcontext,
0534             calibcontext,
0535             &(*pSurface),
0536             ppPlainOptions};
0537
0538     auto trackContainer =
0539         std::make_shared<Acts::VectorTrackContainer>();
0540     auto trackStateContainer =
0541         std::make_shared<Acts::VectorMultiTrajectory>();
0542     ActsTrackFittingAlgorithm::TrackContainer
0543         tracks(trackContainer, trackStateContainer);
0544     auto result = fitTrack(sourceLinks, seed.value(), kfOptions,
0545                            calibrator, tracks);
0546
0547     /// Check that the track fit result did not return an error
0548     if (result.ok())
0549     {
0550       SvtxTrack_v4 newTrack;
0551       newTrack.set_tpc_seed(tpcseed);
0552       newTrack.set_crossing(crossing);
0553       newTrack.set_silicon_seed(siseed);
0554
0555       unsigned int trid = m_trackMap->size();
0556       newTrack.set_id(trid);
0557
0558       if (getTrackFitResult(result, track, &newTrack, tracks, measurements))
0559       {
0560         m_trackMap->insertWithKey(&newTrack, trid);
0561       }
0562       m_nGoodFits++;
0563     }
0564     else
0565     {
0566       m_nBadFits++;
0567       if (Verbosity() > 1)
0568       {
0569         std::cout << "Track fit failed for track " << m_seedMap->find(track)
0570                   << " with Acts error message "
0571                   << result.error() << ", " << result.error().message()
0572                   << std::endl;
0573       }
0574     }
0575   }
0576
0577   return;
0578 }
0579
0580 bool PHCosmicsTrkFitter::getTrackFitResult(FitResult& fitOutput,
0581                                            TrackSeed* seed, SvtxTrack* track,
0582                                            ActsTrackFittingAlgorithm::TrackContainer& tracks,
0583                                            const ActsTrackFittingAlgorithm::MeasurementContainer& measurements)
0584 {
0585   /// Make a trajectory state for storage, which conforms to Acts track fit
0586   /// analysis tool
0587   auto& outtrack = fitOutput.value();
0588   std::vector<Acts::MultiTrajectoryTraits::IndexType> trackTips;
0589   trackTips.reserve(1);
0590   trackTips.emplace_back(outtrack.tipIndex());
0591   Trajectory::IndexedParameters indexedParams;
0592
0593   indexedParams.emplace(std::pair{outtrack.tipIndex(),
0594                                   ActsExamples::TrackParameters{outtrack.referenceSurface().getSharedPtr(),
0595                                                                 outtrack.parameters(), outtrack.covariance(), outtrack.particleHypothesis()}});
0596
0597   if (Verbosity() > 2)
0598   {
0599     std::cout << "Fitted parameters for track" << std::endl;
0600     std::cout << " position : " << outtrack.referenceSurface().localToGlobal(m_transient_geocontext, Acts::Vector2(outtrack.loc0(), outtrack.loc1()), Acts::Vector3(1, 1, 1)).transpose()
0601
0602               << std::endl;
0603
0604     int otcharge = outtrack.qOverP() > 0 ? 1 : -1;
0605     std::cout << "charge: " << otcharge << std::endl;
0606     std::cout << " momentum : " << outtrack.momentum().transpose()
0607               << std::endl;
0608     std::cout << "For trackTip == " << outtrack.tipIndex() << std::endl;
0609   }
0610
0611   Trajectory trajectory(tracks.trackStateContainer(),
0612                         trackTips, indexedParams);
0613
0614   m_trajectories->insert(std::make_pair(track->get_id(), trajectory));
0615
0616   /// Get position, momentum from the Acts output. Update the values of
0617   /// the proto track
0618   updateSvtxTrack(trackTips, indexedParams, tracks, track);
0619
0620   if (m_commissioning)
0621   {
0622     if (track->get_silicon_seed() && track->get_tpc_seed())
0623     {
0624       m_alignStates.fillAlignmentStateMap(tracks, trackTips,
0625                                           track, measurements);
0626     }
0627   }
0628
0629   if (m_actsEvaluator)
0630   {
0631     m_evaluator->evaluateTrackFit(tracks, trackTips, indexedParams, track,
0632                                   seed, measurements);
0633   }
0634
0635   return true;
0636 }
0637
0638 inline ActsTrackFittingAlgorithm::TrackFitterResult PHCosmicsTrkFitter::fitTrack(
0639     const std::vector<Acts::SourceLink>& sourceLinks,
0640     const ActsTrackFittingAlgorithm::TrackParameters& seed,
0641     const ActsTrackFittingAlgorithm::GeneralFitterOptions& kfOptions,
0642     const CalibratorAdapter& calibrator,
0643     ActsTrackFittingAlgorithm::TrackContainer& tracks)
0644 {
0645   return (*m_fitCfg.fit)(sourceLinks, seed, kfOptions, calibrator, tracks);
0646 }
0647
0648 void PHCosmicsTrkFitter::updateSvtxTrack(
0649     std::vector<Acts::MultiTrajectoryTraits::IndexType>& tips,
0650     Trajectory::IndexedParameters& paramsMap,
0651     ActsTrackFittingAlgorithm::TrackContainer& tracks,
0652     SvtxTrack* track)
0653 {
0654   const auto& mj = tracks.trackStateContainer();
0655
0656   /// only one track tip in the track fit Trajectory
0657   auto& trackTip = tips.front();
0658
0659   if (Verbosity() > 2)
0660   {
0661     std::cout << "Identify (proto) track before updating with acts results " << std::endl;
0662     track->identify();
0663   }
0664
0665   // create a state at pathlength = 0.0
0666   // This state holds the track parameters, which will be updated below
0667   float pathlength = 0.0;
0668   SvtxTrackState_v1 out(pathlength);
0669   out.set_x(0.0);
0670   out.set_y(0.0);
0671   out.set_z(0.0);
0672   track->insert_state(&out);
0673
0674   auto trajState =
0675       Acts::MultiTrajectoryHelpers::trajectoryState(mj, trackTip);
0676
0677   const auto& params = paramsMap.find(trackTip)->second;
0678
0679   /// Acts default unit is mm. So convert to cm
0680   track->set_x(params.position(m_transient_geocontext)(0) / Acts::UnitConstants::cm);
0681   track->set_y(params.position(m_transient_geocontext)(1) / Acts::UnitConstants::cm);
0682   track->set_z(params.position(m_transient_geocontext)(2) / Acts::UnitConstants::cm);
0683
0684   track->set_px(params.momentum()(0));
0685   track->set_py(params.momentum()(1));
0686   track->set_pz(params.momentum()(2));
0687
0688   track->set_charge(params.charge());
0689   track->set_chisq(trajState.chi2Sum);
0690   track->set_ndf(trajState.NDF);
0691
0692   ActsTransformations rotater;
0693   rotater.setVerbosity(Verbosity());
0694
0695   if (params.covariance())
0696   {
0697     auto rotatedCov = rotater.rotateActsCovToSvtxTrack(params);
0698
0699     for (int i = 0; i < 6; i++)
0700     {
0701       for (int j = 0; j < 6; j++)
0702       {
0703         track->set_error(i, j, rotatedCov(i, j));
0704       }
0705     }
0706   }
0707
0708   // Also need to update the state list and cluster ID list for all measurements associated with the acts track
0709   // loop over acts track states, copy over to SvtxTrackStates, and add to SvtxTrack
0710
0711   if (m_fillSvtxTrackStates)
0712   {
0713     rotater.fillSvtxTrackStates(mj, trackTip, track,
0714                                 m_transient_geocontext);
0715   }
0716
0717   if (Verbosity() > 2)
0718   {
0719     std::cout << " Identify fitted track after updating track states:"
0720               << std::endl;
0721     track->identify();
0722   }
0723
0724   return;
0725 }
0726
0727 Acts::BoundSquareMatrix PHCosmicsTrkFitter::setDefaultCovariance() const
0728 {
0729   Acts::BoundSquareMatrix cov = Acts::BoundSquareMatrix::Zero();
0730
0731   /// Acts cares about the track covariance as it helps the KF
0732   /// know whether or not to trust the initial track seed or not.
0733   /// We reset it here to some loose values as it helps Acts improve
0734   /// the fitting.
0735   /// If the covariance is too loose, it won't be able to propagate,
0736   /// but if it is too tight, it will just "believe" the track seed over
0737   /// the hit data
0738
0739   // cppcheck-suppress duplicateAssignExpression
0740   double sigmaD0 = 300 * Acts::UnitConstants::um;
0741   double sigmaZ0 = 300 * Acts::UnitConstants::um;
0742   // cppcheck-suppress duplicateAssignExpression
0743   double sigmaPhi = 1 * Acts::UnitConstants::degree;
0744   double sigmaTheta = 1 * Acts::UnitConstants::degree;
0745   double sigmaT = 1. * Acts::UnitConstants::ns;
0746
0747   cov(Acts::eBoundLoc0, Acts::eBoundLoc0) = sigmaD0 * sigmaD0;
0748   cov(Acts::eBoundLoc1, Acts::eBoundLoc1) = sigmaZ0 * sigmaZ0;
0749   cov(Acts::eBoundTime, Acts::eBoundTime) = sigmaT * sigmaT;
0750   cov(Acts::eBoundPhi, Acts::eBoundPhi) = sigmaPhi * sigmaPhi;
0751   cov(Acts::eBoundTheta, Acts::eBoundTheta) = sigmaTheta * sigmaTheta;
0752   /// Acts takes this value very seriously - tuned to be in a "sweet spot"
0753   cov(Acts::eBoundQOverP, Acts::eBoundQOverP) = 0.0001;
0754
0755   return cov;
0756 }
0757
0758 void PHCosmicsTrkFitter::printTrackSeed(const ActsTrackFittingAlgorithm::TrackParameters& seed) const
0759 {
0760   std::cout
0761       << PHWHERE
0762       << " Processing proto track:"
0763       << std::endl;
0764
0765   std::cout
0766       << "position: " << seed.position(m_transient_geocontext).transpose()
0767       << std::endl
0768       << "momentum: " << seed.momentum().transpose()
0769       << std::endl;
0770
0771   std::cout << "charge : " << seed.charge() << std::endl;
0772   std::cout << "absolutemom : " << seed.absoluteMomentum() << std::endl;
0773 }
0774
0775 int PHCosmicsTrkFitter::createNodes(PHCompositeNode* topNode)
0776 {
0777   PHNodeIterator iter(topNode);
0778
0779   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0780
0781   if (!dstNode)
0782   {
0783     std::cerr << "DST node is missing, quitting" << std::endl;
0784     throw std::runtime_error("Failed to find DST node in PHCosmicsTrkFitter::createNodes");
0785   }
0786
0787   PHNodeIterator dstIter(topNode);
0788   PHCompositeNode* svtxNode = dynamic_cast<PHCompositeNode*>(dstIter.findFirst("PHCompositeNode", "SVTX"));
0789
0790   if (!svtxNode)
0791   {
0792     svtxNode = new PHCompositeNode("SVTX");
0793     dstNode->addNode(svtxNode);
0794   }
0795
0796   m_trajectories = findNode::getClass<std::map<const unsigned int, Trajectory>>(topNode, "ActsTrajectories");
0797   if (!m_trajectories)
0798   {
0799     m_trajectories = new std::map<const unsigned int, Trajectory>;
0800     auto node =
0801         new PHDataNode<std::map<const unsigned int, Trajectory>>(m_trajectories, "ActsTrajectories");
0802     svtxNode->addNode(node);
0803   }
0804
0805   m_trackMap = findNode::getClass<SvtxTrackMap>(topNode, _track_map_name);
0806
0807   if (!m_trackMap)
0808   {
0809     m_trackMap = new SvtxTrackMap_v2;
0810     PHIODataNode<PHObject>* node = new PHIODataNode<PHObject>(m_trackMap, _track_map_name, "PHObject");
0811     svtxNode->addNode(node);
0812   }
0813
0814   m_alignmentStateMap = findNode::getClass<SvtxAlignmentStateMap>(topNode, "SvtxAlignmentStateMap");
0815   if (!m_alignmentStateMap)
0816   {
0817     m_alignmentStateMap = new SvtxAlignmentStateMap_v1;
0818     auto node = new PHDataNode<SvtxAlignmentStateMap>(m_alignmentStateMap, "SvtxAlignmentStateMap", "PHObject");
0819     svtxNode->addNode(node);
0820   }
0821
0822   m_alignmentTransformationMapTransient = findNode::getClass<alignmentTransformationContainer>(topNode, "alignmentTransformationContainerTransient");
0823   if (!m_alignmentTransformationMapTransient)
0824   {
0825     std::cout << PHWHERE << "alignmentTransformationContainerTransient not on node tree. Bailing"
0826               << std::endl;
0827     return Fun4AllReturnCodes::ABORTEVENT;
0828   }
0829
0830   if (m_actsEvaluator)
0831   {
0832     m_seedTracks = findNode::getClass<SvtxTrackMap>(topNode, _seed_track_map_name);
0833
0834     if (!m_seedTracks)
0835     {
0836       m_seedTracks = new SvtxTrackMap_v2;
0837
0838       PHIODataNode<PHObject>* seedNode =
0839           new PHIODataNode<PHObject>(m_seedTracks, _seed_track_map_name, "PHObject");
0840       svtxNode->addNode(seedNode);
0841     }
0842   }
0843
0844   return Fun4AllReturnCodes::EVENT_OK;
0845 }
0846
0847 int PHCosmicsTrkFitter::getNodes(PHCompositeNode* topNode)
0848 {
0849   // tpc seeds
0850   m_tpcSeeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
0851   if (!m_tpcSeeds)
0852   {
0853     std::cout << PHWHERE << "TpcTrackSeedContainer not on node tree. Bailing"
0854               << std::endl;
0855     return Fun4AllReturnCodes::ABORTEVENT;
0856   }
0857
0858   // silicon seeds
0859   m_siliconSeeds = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
0860   if (!m_siliconSeeds)
0861   {
0862     std::cout << PHWHERE << "SiliconTrackSeedContainer not on node tree. Bailing"
0863               << std::endl;
0864     return Fun4AllReturnCodes::ABORTEVENT;
0865   }
0866
0867   // clusters
0868   m_clusterContainer = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0869   if (!m_clusterContainer)
0870   {
0871     std::cout << PHWHERE
0872               << "No trkr cluster container, exiting." << std::endl;
0873     return Fun4AllReturnCodes::ABORTEVENT;
0874   }
0875
0876   // acts geometry
0877   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0878   if (!m_tGeometry)
0879   {
0880     std::cout << "ActsGeometry not on node tree. Exiting."
0881               << std::endl;
0882
0883     return Fun4AllReturnCodes::ABORTEVENT;
0884   }
0885
0886   // tracks
0887   m_seedMap = findNode::getClass<TrackSeedContainer>(topNode, "SvtxTrackSeedContainer");
0888   if (!m_seedMap)
0889   {
0890     std::cout << "No Svtx seed map on node tree. Exiting."
0891               << std::endl;
0892     return Fun4AllReturnCodes::ABORTEVENT;
0893   }
0894
0895   // tpc global position wrapper
0896   m_globalPositionWrapper.loadNodes(topNode);
0897
0898   return Fun4AllReturnCodes::EVENT_OK;
0899 }
0900
0901 void PHCosmicsTrkFitter::makeBranches()
0902 {
0903   m_tree->Branch("seed", &m_seed, "m_seed/I");
0904   m_tree->Branch("event", &m_event, "m_event/I");
0905   m_tree->Branch("R", &m_R, "m_R/F");
0906   m_tree->Branch("X0", &m_X0, "m_X0/F");
0907   m_tree->Branch("Y0", &m_Y0, "m_Y0/F");
0908   m_tree->Branch("Z0", &m_Z0, "m_Z0/F");
0909   m_tree->Branch("slope", &m_slope, "m_slope/F");
0910   m_tree->Branch("pcax", &m_pcax, "m_pcax/F");
0911   m_tree->Branch("pcay", &m_pcay, "m_pcay/F");
0912   m_tree->Branch("pcaz", &m_pcaz, "m_pcaz/F");
0913   m_tree->Branch("px", &m_px, "m_px/F");
0914   m_tree->Branch("py", &m_py, "m_py/F");
0915   m_tree->Branch("pz", &m_pz, "m_pz/F");
0916   m_tree->Branch("charge", &m_charge, "m_charge/I");
0917   m_tree->Branch("nmaps", &m_nmaps, "m_nmaps/I");
0918   m_tree->Branch("nintt", &m_nintt, "m_nintt/I");
0919   m_tree->Branch("ntpc", &m_ntpc, "m_ntpc/I");
0920   m_tree->Branch("nmm", &m_nmm, "m_nmm/I");
0921   m_tree->Branch("locx", &m_locx);
0922   m_tree->Branch("locy", &m_locy);
0923   m_tree->Branch("x", &m_x);
0924   m_tree->Branch("y", &m_y);
0925   m_tree->Branch("z", &m_z);
0926   m_tree->Branch("r", &m_r);
0927   m_tree->Branch("layer", &m_layer);
0928   m_tree->Branch("phi", &m_phi);
0929   m_tree->Branch("eta", &m_eta);
0930   m_tree->Branch("phisize", &m_phisize);
0931   m_tree->Branch("zsize", &m_zsize);
0932   m_tree->Branch("ephi", &m_ephi);
0933   m_tree->Branch("ez", &m_ez);
0934 }
0935 void PHCosmicsTrkFitter::fillVectors(TrackSeed* tpcseed, TrackSeed* siseed)
0936 {
0937   for (auto seed : {tpcseed, siseed})
0938   {
0939     if (!seed)
0940     {
0941       continue;
0942     }
0943
0944     for (auto it = seed->begin_cluster_keys(); it != seed->end_cluster_keys();
0945          ++it)
0946     {
0947       auto key = *it;
0948       auto cluster = m_clusterContainer->findCluster(key);
0949       m_locx.push_back(cluster->getLocalX());
0950       float ly = cluster->getLocalY();
0951       if (TrkrDefs::getTrkrId(key) == TrkrDefs::TrkrId::tpcId)
0952       {
0953         double drift_velocity = m_tGeometry->get_drift_velocity();
0954         double zdriftlength = cluster->getLocalY() * drift_velocity;
0955         double surfCenterZ = 52.89;                // 52.89 is where G4 thinks the surface center is
0956         double zloc = surfCenterZ - zdriftlength;  // converts z drift length to local z position in the TPC in north
0957         unsigned int side = TpcDefs::getSide(key);
0958         if (side == 0)
0959         {
0960           zloc = -zloc;
0961         }
0962         ly = zloc * 10;
0963       }
0964       m_locy.push_back(ly);
0965       auto glob = m_tGeometry->getGlobalPosition(key, cluster);
0966       m_x.push_back(glob.x());
0967       m_y.push_back(glob.y());
0968       m_z.push_back(glob.z());
0969       float r = std::sqrt(glob.x() * glob.x() + glob.y() * glob.y());
0970       m_r.push_back(r);
0971       TVector3 globt(glob.x(), glob.y(), glob.z());
0972       m_phi.push_back(globt.Phi());
0973       m_eta.push_back(globt.Eta());
0974       m_phisize.push_back(cluster->getPhiSize());
0975       m_zsize.push_back(cluster->getZSize());
0976       auto para_errors =
0977           m_clusErrPara.get_clusterv5_modified_error(cluster, r, key);
0978
0979       m_ephi.push_back(std::sqrt(para_errors.first));
0980       m_ez.push_back(std::sqrt(para_errors.second));
0981     }
0982   }
0983 }
0984 void PHCosmicsTrkFitter::clearVectors()
0985 {
0986   m_locx.clear();
0987   m_locy.clear();
0988   m_x.clear();
0989   m_y.clear();
0990   m_z.clear();
0991   m_layer.clear();
0992   m_r.clear();
0993   m_phi.clear();
0994   m_eta.clear();
0995   m_phisize.clear();
0996   m_zsize.clear();
0997   m_ephi.clear();
0998   m_ez.clear();
0999 }
1000
1001 void PHCosmicsTrkFitter::getCharge(
1002     TrackSeed* track,
1003     // TrkrClusterContainer*  clusterContainer,
1004     // ActsGeometry* tGeometry,
1005     // alignmentTransformationContainer* transformMapTransient,
1006     // float vertexRadius,
1007     int& charge,
1008     float& cosmicslope)
1009 {
1010   Acts::GeometryContext transient_geocontext;
1011   transient_geocontext = m_alignmentTransformationMapTransient;  // set local/global transforms to distortion corrected ones for this track
1012
1013   std::vector<Acts::Vector3> global_vec;
1014
1015   for (auto clusIter = track->begin_cluster_keys();
1016        clusIter != track->end_cluster_keys();
1017        ++clusIter)
1018   {
1019     auto key = *clusIter;
1020     auto cluster = m_clusterContainer->findCluster(key);
1021     if (!cluster)
1022     {
1023       std::cout << "MakeSourceLinks::getCharge: Failed to get cluster with key " << key << " for track seed" << std::endl;
1024       continue;
1025     }
1026
1027     auto surf = m_tGeometry->maps().getSurface(key, cluster);
1028     if (!surf)
1029     {
1030       continue;
1031     }
1032
1033     // get cluster global positions
1034     Acts::Vector2 local = m_tGeometry->getLocalCoords(key, cluster);  // converts TPC time to z
1035     Acts::Vector3 glob = surf->localToGlobal(transient_geocontext,
1036                                              local * Acts::UnitConstants::cm,
1037                                              Acts::Vector3(1, 1, 1));
1038     glob /= Acts::UnitConstants::cm;
1039
1040     global_vec.push_back(glob);
1041   }
1042
1043   Acts::Vector3 globalMostOuter(std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN(), std::numeric_limits<double>::quiet_NaN());
1044   Acts::Vector3 globalSecondMostOuter(0, 999999, 0);
1045   float largestR = 0;
1046   // loop over global positions
1047   for (auto& i : global_vec)
1048   {
1049     Acts::Vector3 global = i;
1050     // float r = std::sqrt(square(global.x()) + square(global.y()));
1051     float r = radius(global.x(), global.y());
1052
1053     /// use the top hemisphere to determine the charge
1054     if (r > largestR && global.y() > 0)
1055     {
1056       globalMostOuter = i;
1057       largestR = r;
1058     }
1059   }
1060
1061   //! find the closest cluster to the outermost cluster
1062   float maxdr = std::numeric_limits<float>::max();
1063   for (auto& i : global_vec)
1064   {
1065     if (i.y() < 0)
1066     {
1067       continue;
1068     }
1069
1070     float dr = std::sqrt(square(globalMostOuter.x()) + square(globalMostOuter.y())) - std::sqrt(square(i.x()) + square(i.y()));
1071     //! Place a dr cut to get maximum bend due to TPC clusters having
1072     //! larger fluctuations
1073     if (dr < maxdr && dr > 10)
1074     {
1075       maxdr = dr;
1076       globalSecondMostOuter = i;
1077     }
1078   }
1079
1080   //! we have to calculate phi WRT the vertex position outside the detector,
1081   //! not at (0,0)
1082   Acts::Vector3 vertex(0, m_vertexRadius, 0);
1083   globalMostOuter -= vertex;
1084   globalSecondMostOuter -= vertex;
1085
1086   const auto firstphi = atan2(globalMostOuter.y(), globalMostOuter.x());
1087   const auto secondphi = atan2(globalSecondMostOuter.y(),
1088                                globalSecondMostOuter.x());
1089   auto dphi = secondphi - firstphi;
1090
1091   if (dphi > M_PI)
1092   {
1093     dphi = 2. * M_PI - dphi;
1094   }
1095   if (dphi < -M_PI)
1096   {
1097     dphi = 2 * M_PI + dphi;
1098   }
1099
1100   if (dphi > 0)
1101   {
1102     charge = -1;
1103   }
1104   else
1105   {
1106     charge = 1;
1107   }
1108
1109   float r1 = std::sqrt(square(globalMostOuter.x()) + square(globalMostOuter.y()));
1110   float r2 = std::sqrt(square(globalSecondMostOuter.x()) + square(globalSecondMostOuter.y()));
1111   float z1 = globalMostOuter.z();
1112   float z2 = globalSecondMostOuter.z();
1113   cosmicslope = (r2 - r1) / (z2 - z1);
1114
1115   return;
1116 }