From 9d208c1de8f55061c169852eff5bdc4ec6b1f1e9 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Fri, 10 Jun 2022 17:22:41 +0200 Subject: [PATCH 01/14] add 'magnetic_field' parameter in Track and Tracking4D --- src/modules/Tracking4D/Tracking4D.cpp | 3 +++ src/modules/Tracking4D/Tracking4D.h | 2 ++ src/objects/Track.cpp | 4 ++++ src/objects/Track.hpp | 7 +++++++ 4 files changed, 16 insertions(+) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index 5903f0d4..1fe15938 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -34,6 +34,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("volume_scattering", false); config_.setDefault("reject_by_roi", false); config_.setDefault("unique_cluster_usage", false); + config_.setDefault("magnetic_field", ROOT::Math::XYZVector()); if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { config_.setDefault("time_cut_rel", 3.0); @@ -66,6 +67,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("volume_scattering"); reject_by_ROI_ = config_.get("reject_by_roi"); unique_cluster_usage_ = config_.get("unique_cluster_usage"); + magnetic_field_ = config_.get("magnetic_field"); // print a warning if volumeScatterer are used as this causes fit failures // that are still not understood @@ -322,6 +324,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { track->setVolumeScatter(volume_radiation_length_); } track->setParticleMomentum(momentum_); + track->setMagneticField(magnetic_field_); // Loop over each subsequent plane and look for a cluster within the timing cuts size_t detector_nr = 2; diff --git a/src/modules/Tracking4D/Tracking4D.h b/src/modules/Tracking4D/Tracking4D.h index 91ad1acd..6415c70b 100644 --- a/src/modules/Tracking4D/Tracking4D.h +++ b/src/modules/Tracking4D/Tracking4D.h @@ -14,6 +14,7 @@ #include #include #include +#include #include #include "core/module/Module.hpp" #include "objects/Cluster.hpp" @@ -92,6 +93,7 @@ namespace corryvreckan { std::map, XYVector> spatial_cuts_; std::string timestamp_from_; std::string track_model_; + ROOT::Math::XYZVector magnetic_field_; // Function to calculate the weighted average timestamp from the clusters of a track double calculate_average_timestamp(const Track* track); diff --git a/src/objects/Track.cpp b/src/objects/Track.cpp index d21c07a0..6449c4d8 100644 --- a/src/objects/Track.cpp +++ b/src/objects/Track.cpp @@ -124,6 +124,10 @@ void Track::setParticleMomentum(double p) { momentum_ = p; } +void Track::setMagneticField(ROOT::Math::XYZVector B) { + magneticField_ = B; +} + double Track::getChi2() const { if(!isFitted_) { throw RequestParameterBeforeFitError(this, "chi2"); diff --git a/src/objects/Track.hpp b/src/objects/Track.hpp index c1bbdfa2..549dbefb 100644 --- a/src/objects/Track.hpp +++ b/src/objects/Track.hpp @@ -105,6 +105,12 @@ namespace corryvreckan { */ void setParticleMomentum(double momentum); + /** + * @brief Set the magnetic field + * @param magneticField + */ + void setMagneticField(ROOT::Math::XYZVector magneticField); + /** * @brief Get the chi2 of the track fit * @return chi2 @@ -306,6 +312,7 @@ namespace corryvreckan { double chi2ndof_; bool isFitted_{}; double momentum_{-1}; + ROOT::Math::XYZVector magneticField_; // ROOT I/O class definition - update version number when you change this class! ClassDefOverride(Track, 11) -- GitLab From eec4f1596c5d4dcfc137a4db0d561d6bdb1f4bc3 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Fri, 17 Jun 2022 13:54:09 +0200 Subject: [PATCH 02/14] Define Standard units to be compatible with the manual --- src/core/Corryvreckan.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/core/Corryvreckan.cpp b/src/core/Corryvreckan.cpp index 50705e63..61817f97 100644 --- a/src/core/Corryvreckan.cpp +++ b/src/core/Corryvreckan.cpp @@ -238,10 +238,11 @@ void Corryvreckan::add_units() { // NOTE: fixed by above Units::add("V", 1e-6); Units::add("kV", 1e-3); + Units::add("MV", 1); // MAGNETIC FIELD - Units::add("T", 1e-3); - Units::add("mT", 1e-6); + Units::add("T", 1); + Units::add("mT", 1e-3); // ANGLES // NOTE: these are fake units -- GitLab From a291c9702b7088c46a76850d4e7bad330012a338 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Fri, 17 Jun 2022 15:07:21 +0200 Subject: [PATCH 03/14] Store all coordinates of origin in Track::Plane instead of z only --- .../AlignmentTrackChi2/AlignmentTrackChi2.cpp | 2 +- .../AnalysisTelescope/AnalysisTelescope.cpp | 80 ++++++++++++------- .../AnalysisTelescope/AnalysisTelescope.h | 12 ++- src/modules/Tracking4D/Tracking4D.cpp | 2 +- .../TrackingMultiplet/TrackingMultiplet.cpp | 2 +- src/objects/GblTrack.cpp | 10 +-- src/objects/Track.cpp | 16 ++-- src/objects/Track.hpp | 11 +-- 8 files changed, 82 insertions(+), 53 deletions(-) diff --git a/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp b/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp index 5e752386..21579c77 100644 --- a/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp +++ b/src/modules/AlignmentTrackChi2/AlignmentTrackChi2.cpp @@ -133,7 +133,7 @@ void AlignmentTrackChi2::MinimiseTrackChi2(Int_t&, Double_t*, Double_t& result, // Refit the track track->registerPlane(AlignmentTrackChi2::globalDetector->getName(), - AlignmentTrackChi2::globalDetector->displacement().z(), + AlignmentTrackChi2::globalDetector->origin(), AlignmentTrackChi2::globalDetector->materialBudget(), AlignmentTrackChi2::globalDetector->toLocal()); LOG(DEBUG) << "Updated transformations for detector " << AlignmentTrackChi2::globalDetector->getName(); diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp index 67b50d38..e388ce87 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.cpp +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.cpp @@ -38,27 +38,35 @@ void AnalysisTelescope::initialize() { if(detector->isDUT()) { std::string title = detector->getName() + " Telescope resolution X;x_{track}-x_{MC} [mm];events"; - telescopeResolutionX[detector->getName()] = new TH1F("telescopeResolutionX", title.c_str(), 600, -0.2, 0.2); + telescopeResolutionX[detector->getName()] = new TH1F("telescopeResolutionX", title.c_str(), 600, -3*detector->getPitch().X(), 3*detector->getPitch().X()); title = detector->getName() + " Telescope resolution Y;y_{track}-y_{MC} [mm];events"; - telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -0.2, 0.2); + telescopeResolutionY[detector->getName()] = new TH1F("telescopeResolutionY", title.c_str(), 600, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); } else { std::string title = detector->getName() + " Biased residual X (local);x_{track}-x_{cluster} [mm];events"; - telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -0.2, 0.2); + telescopeResidualsLocalX[detector->getName()] = new TH1F("residualX_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); title = detector->getName() + " Biased residual Y (local);y_{track}-y_{cluster} [mm];events"; - telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -0.2, 0.2); + telescopeResidualsLocalY[detector->getName()] = new TH1F("residualY_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); title = detector->getName() + " Biased residual X (global);x_{track}-x_{cluster} [mm];events"; telescopeResidualsX[detector->getName()] = new TH1F("residualX_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " Biased residual Y (global);y_{track}-y_{cluster} [mm];events"; telescopeResidualsY[detector->getName()] = new TH1F("residualY_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual X (local);x_{track}-x_{MC} [mm];events"; - telescopeMCresidualsLocalX[detector->getName()] = new TH1F("residualX_MC_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual Y (local);y_{track}-y_{MC} [mm];events"; - telescopeMCresidualsLocalY[detector->getName()] = new TH1F("residualY_MC_local", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual X (global);x_{track}-x_{MC} [mm];events"; - telescopeMCresidualsX[detector->getName()] = new TH1F("residualX_MC_global", title.c_str(), 400, -0.2, 0.2); - title = detector->getName() + " Biased MC residual Y (global);y_{track}-y_{MC} [mm];events"; - telescopeMCresidualsY[detector->getName()] = new TH1F("residualY_MC_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC cluster residual X (local);x_{cluster}-x_{MC} [mm];events"; + telescopeMCresidualsClusterLocalX[detector->getName()] = new TH1F("residualX_MC_cluster_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); + title = detector->getName() + " Biased MC cluster residual Y (local);y_{cluster}-y_{MC} [mm];events"; + telescopeMCresidualsClusterLocalY[detector->getName()] = new TH1F("residualY_MC_cluster_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); + title = detector->getName() + " Biased MC track residual X (local);x_{track}-x_{MC} [mm];events"; + telescopeMCresidualsTrackLocalX[detector->getName()] = new TH1F("residualX_MC_track_local", title.c_str(), 400, -3*detector->getPitch().X(), 3*detector->getPitch().X()); + title = detector->getName() + " Biased MC track residual Y (local);y_{track}-y_{MC} [mm];events"; + telescopeMCresidualsTrackLocalY[detector->getName()] = new TH1F("residualY_MC_track_local", title.c_str(), 400, -3*detector->getPitch().Y(), 3*detector->getPitch().Y()); + title = detector->getName() + " Biased MC cluster residual X (global);x_{cluster}-x_{MC} [mm];events"; + telescopeMCresidualsClusterX[detector->getName()] = new TH1F("residualX_MC_cluster_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC cluster residual Y (global);y_{cluster}-y_{MC} [mm];events"; + telescopeMCresidualsClusterY[detector->getName()] = new TH1F("residualY_MC_cluster_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC track residual X (global);x_{track}-x_{MC} [mm];events"; + telescopeMCresidualsTrackX[detector->getName()] = new TH1F("residualX_MC_track_global", title.c_str(), 400, -0.2, 0.2); + title = detector->getName() + " Biased MC track residual Y (global);y_{track}-y_{MC} [mm];events"; + telescopeMCresidualsTrackY[detector->getName()] = new TH1F("residualY_MC_track_global", title.c_str(), 400, -0.2, 0.2); title = detector->getName() + " pixel - seed timestamp (all pixels w/o seed);ts_{pixel} - ts_{seed} [ns]; pixel charge [e];events"; @@ -145,7 +153,7 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { if(track->getType() != "GblTrack") { intercept = track->getIntercept(cluster->global().z()); } - auto interceptLocal = detector->globalToLocal(intercept); + auto interceptLocal = detector->getLocalIntercept(track.get()); telescopeResidualsLocalX[name]->Fill(cluster->local().x() - interceptLocal.X()); telescopeResidualsLocalY[name]->Fill(cluster->local().y() - interceptLocal.Y()); telescopeResidualsX[name]->Fill(cluster->global().x() - intercept.X()); @@ -179,13 +187,24 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { continue; } - ROOT::Math::XYZPoint particlePosition = closestApproach(cluster->local(), mcParticles); - telescopeMCresidualsLocalX[name]->Fill(cluster->local().x() + detector->getSize().X() / 2 - - particlePosition.X()); - telescopeMCresidualsLocalY[name]->Fill(cluster->local().y() + detector->getSize().Y() / 2 - - particlePosition.Y()); - telescopeMCresidualsX[name]->Fill(interceptLocal.X() + detector->getSize().X() / 2 - particlePosition.X()); - telescopeMCresidualsY[name]->Fill(interceptLocal.Y() + detector->getSize().Y() / 2 - particlePosition.Y()); + ROOT::Math::XYZPoint clusterPosLocalConv( + cluster->local().X()+detector->getSize().X()/2-detector->getPitch().X()/2, + cluster->local().Y()+detector->getSize().Y()/2-detector->getPitch().Y()/2, + cluster->local().Z()); + ROOT::Math::XYZPoint particlePositionConv = closestApproach(clusterPosLocalConv, mcParticles); + ROOT::Math::XYZPoint particlePosition( + particlePositionConv.X() - detector->getSize().X()/2 + detector->getPitch().X()/2, + particlePositionConv.Y() - detector->getSize().Y()/2 + detector->getPitch().Y()/2, + particlePositionConv.Z()); + ROOT::Math::XYZPoint particlePositionGlobal = detector->localToGlobal(particlePosition); + telescopeMCresidualsClusterLocalX[name]->Fill(cluster->local().X() - particlePosition.X()); + telescopeMCresidualsClusterLocalY[name]->Fill(cluster->local().Y() - particlePosition.Y()); + telescopeMCresidualsTrackLocalX[name]->Fill(interceptLocal.X()- particlePosition.X()); + telescopeMCresidualsTrackLocalY[name]->Fill(interceptLocal.Y()- particlePosition.Y()); + telescopeMCresidualsClusterX[name]->Fill(cluster->global().X() - particlePositionGlobal.X()); + telescopeMCresidualsClusterY[name]->Fill(cluster->global().X() - particlePositionGlobal.X()); + telescopeMCresidualsTrackX[name]->Fill(intercept.X() - particlePositionGlobal.X()); + telescopeMCresidualsTrackY[name]->Fill(intercept.X() - particlePositionGlobal.X()); } // Calculate telescope resolution at DUTs @@ -196,14 +215,19 @@ StatusCode AnalysisTelescope::run(const std::shared_ptr& clipboard) { continue; } - auto intercept = detector->getIntercept(track.get()); - auto interceptLocal = detector->globalToLocal(intercept); - auto particlePosition = closestApproach(interceptLocal, mcParticles); - - telescopeResolutionX[detector->getName()]->Fill(interceptLocal.X() + detector->getSize().X() / 2 - - particlePosition.X()); - telescopeResolutionY[detector->getName()]->Fill(interceptLocal.Y() + detector->getSize().Y() / 2 - - particlePosition.Y()); + auto interceptLocal = detector->getLocalIntercept(track.get()); + ROOT::Math::XYZPoint interceptLocalConv( + interceptLocal.X()+detector->getSize().X()/2-detector->getPitch().X()/2, + interceptLocal.Y()+detector->getSize().Y()/2-detector->getPitch().Y()/2, + interceptLocal.Z()); + auto particlePositionConv = closestApproach(interceptLocalConv, mcParticles); + ROOT::Math::XYZPoint particlePosition( + particlePositionConv.X() - detector->getSize().X()/2 + detector->getPitch().X()/2, + particlePositionConv.Y() - detector->getSize().Y()/2 + detector->getPitch().Y()/2, + particlePositionConv.Z()); + + telescopeResolutionX[detector->getName()]->Fill(interceptLocal.X() - particlePosition.X()); + telescopeResolutionY[detector->getName()]->Fill(interceptLocal.Y() - particlePosition.Y()); } } diff --git a/src/modules/AnalysisTelescope/AnalysisTelescope.h b/src/modules/AnalysisTelescope/AnalysisTelescope.h index 30e52fa4..671ac86a 100644 --- a/src/modules/AnalysisTelescope/AnalysisTelescope.h +++ b/src/modules/AnalysisTelescope/AnalysisTelescope.h @@ -35,10 +35,14 @@ namespace corryvreckan { ROOT::Math::XYZPoint closestApproach(ROOT::Math::XYZPoint position, const MCParticleVector& particles); // Histograms for each of the devices - std::map telescopeMCresidualsLocalX; - std::map telescopeMCresidualsLocalY; - std::map telescopeMCresidualsX; - std::map telescopeMCresidualsY; + std::map telescopeMCresidualsClusterLocalX; + std::map telescopeMCresidualsClusterLocalY; + std::map telescopeMCresidualsTrackLocalX; + std::map telescopeMCresidualsTrackLocalY; + std::map telescopeMCresidualsClusterX; + std::map telescopeMCresidualsClusterY; + std::map telescopeMCresidualsTrackX; + std::map telescopeMCresidualsTrackY; std::map telescopeResidualsLocalX; std::map telescopeResidualsLocalY; diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index 1fe15938..d764ed08 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -336,7 +336,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { auto detectorID = detector->getName(); LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); track->registerPlane( - detectorID, detector->displacement().z(), detector->materialBudget(), detector->toLocal()); + detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); if(detector == reference_first || detector == reference_last) { continue; diff --git a/src/modules/TrackingMultiplet/TrackingMultiplet.cpp b/src/modules/TrackingMultiplet/TrackingMultiplet.cpp index 5fda24f4..1fcbe31a 100644 --- a/src/modules/TrackingMultiplet/TrackingMultiplet.cpp +++ b/src/modules/TrackingMultiplet/TrackingMultiplet.cpp @@ -294,7 +294,7 @@ TrackVector TrackingMultiplet::find_multiplet_tracklets(const streams& stream, for(auto& det : cluster_trees) { auto detector = det.first.get(); trackletCandidate->registerPlane( - detector->getName(), detector->displacement().x(), detector->materialBudget(), detector->toLocal()); + detector->getName(), detector->origin(), detector->materialBudget(), detector->toLocal()); } trackletCandidate->addCluster(clusterFirst.get()); trackletCandidate->addCluster(clusterLast.get()); diff --git a/src/objects/GblTrack.cpp b/src/objects/GblTrack.cpp index b80eb7c0..a37c84e1 100644 --- a/src/objects/GblTrack.cpp +++ b/src/objects/GblTrack.cpp @@ -211,7 +211,7 @@ void GblTrack::prepare_gblpoints() { for(auto& l : planes_) { total_material += l.getMaterialBudget(); if(clusters.count(l.getName()) == 1) { - l.setPosition(clusters.at(l.getName())->global().z()); + l.setPosition(clusters.at(l.getName())->global()); l.setCluster(clusters.at(l.getName())); } } @@ -221,7 +221,7 @@ void GblTrack::prepare_gblpoints() { // add volume scattering length - for now simply the distance between first and last plane if(use_volume_scatter_) { - total_material += (planes_.back().getPosition() - planes_.front().getPosition()) / scattering_length_volume_; + total_material += (planes_.back().getPosition().z() - planes_.front().getPosition().z()) / scattering_length_volume_; } // Prepare transformations for first plane: @@ -355,7 +355,7 @@ ROOT::Math::XYZPoint GblTrack::getIntercept(double z) const { } for(const auto& l : planes_) { - if(l.getPosition() >= z) { + if(l.getPosition().z() >= z) { found = true; break; } @@ -408,7 +408,7 @@ XYZPoint GblTrack::get_position_outside_telescope(double z) const { auto last = planes_.end(); last--; // No direct iterator for last element // check if z is up or downstream - bool upstream = (z < first->getPosition()); + bool upstream = (z < first->getPosition().z()); auto outerPlane = (upstream ? first->getName() : last->getName()); // inner neighbour of plane - simply adjust the iterators @@ -447,7 +447,7 @@ ROOT::Math::XYZVector GblTrack::getDirection(const std::string& detectorID) cons } XYZVector GblTrack::getDirection(const double& z) const { - auto planeUpstream = std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition() > z; }); + auto planeUpstream = std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition().z() > z; }); if(planeUpstream != planes_.end()) { return getDirection(planeUpstream->getName()); } else { diff --git a/src/objects/Track.cpp b/src/objects/Track.cpp index 6449c4d8..11932f1d 100644 --- a/src/objects/Track.cpp +++ b/src/objects/Track.cpp @@ -13,11 +13,11 @@ using namespace corryvreckan; -Track::Plane::Plane(std::string name, double z, double x_x0, Transform3D to_local) - : z_(z), x_x0_(x_x0), name_(std::move(name)), to_local_(to_local) {} +Track::Plane::Plane(std::string name, ROOT::Math::XYZPoint origin, double x_x0, Transform3D to_local) + : origin_(origin), x_x0_(x_x0), name_(std::move(name)), to_local_(to_local) {} -double Track::Plane::getPosition() const { - return z_; +ROOT::Math::XYZPoint Track::Plane::getPosition() const { + return origin_; } double Track::Plane::getMaterialBudget() const { @@ -49,7 +49,7 @@ Transform3D Track::Plane::getToGlobal() const { } bool Track::Plane::operator<(const Plane& pl) const { - return z_ < pl.z_; + return origin_.z() < pl.origin_.z(); } void Track::Plane::setCluster(const Cluster* cluster) { @@ -57,7 +57,7 @@ void Track::Plane::setCluster(const Cluster* cluster) { } void Track::Plane::print(std::ostream& os) const { - os << "Plane at " << z_ << " with rad. length " << x_x0_ << ", name " << name_ << " and"; + os << "Plane at " << origin_ << " with rad. length " << x_x0_ << ", name " << name_ << " and"; if(hasCluster()) { os << "cluster with global pos: " << getCluster()->global(); } else { @@ -229,8 +229,8 @@ double Track::getMaterialBudget(const std::string& detectorID) const { return budget; } -void Track::registerPlane(const std::string& name, double z, double x0, Transform3D g2l) { - Plane p(name, z, x0, g2l); +void Track::registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, double x0, Transform3D g2l) { + Plane p(name, origin, x0, g2l); auto pl = std::find_if(planes_.begin(), planes_.end(), [&p](const Plane& plane) { return plane.getName() == p.getName(); }); if(pl == planes_.end()) { diff --git a/src/objects/Track.hpp b/src/objects/Track.hpp index 549dbefb..9b799ed5 100644 --- a/src/objects/Track.hpp +++ b/src/objects/Track.hpp @@ -239,12 +239,12 @@ namespace corryvreckan { virtual void setVolumeScatter(double length) = 0; - void registerPlane(const std::string& name, double z, double x0, Transform3D g2l); + void registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, double x0, Transform3D g2l); class Plane { public: Plane() = default; - Plane(std::string name, double z, double x_x0, Transform3D to_local); + Plane(std::string name, ROOT::Math::XYZPoint origin, double x_x0, Transform3D to_local); /** * @brief Required virtual destructor */ @@ -266,7 +266,7 @@ namespace corryvreckan { Plane& operator=(Plane&& rhs) = default; /// @} - double getPosition() const; + ROOT::Math::XYZPoint getPosition() const; double getMaterialBudget() const; bool hasCluster() const; @@ -278,7 +278,7 @@ namespace corryvreckan { // sorting overload bool operator<(const Plane& pl) const; // set elements that might be unknown at construction - void setPosition(double z) { z_ = z; } + void setPosition(ROOT::Math::XYZPoint origin) { origin_ = origin; } void setCluster(const Cluster* cluster); void print(std::ostream& os) const; @@ -289,7 +289,8 @@ namespace corryvreckan { void petrifyHistory(); private: - double z_, x_x0_; + ROOT::Math::XYZPoint origin_; + double x_x0_; std::string name_; PointerWrapper cluster_; Transform3D to_local_; -- GitLab From f0f583285b8b5b65304bd50b07a809faac96ed42 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 21 Jun 2022 11:42:27 +0200 Subject: [PATCH 04/14] Add new track model HelixTrack --- src/core/detector/PixelDetector.cpp | 2 +- src/objects/CMakeLists.txt | 2 + src/objects/HelixTrack.cpp | 383 ++++++++++++++++++++++++++++ src/objects/HelixTrack.hpp | 125 +++++++++ src/objects/Linkdef.h | 1 + src/objects/Track.cpp | 4 +- src/objects/Track.hpp | 3 +- src/objects/objects.h | 2 +- 8 files changed, 518 insertions(+), 4 deletions(-) create mode 100644 src/objects/HelixTrack.cpp create mode 100644 src/objects/HelixTrack.hpp diff --git a/src/core/detector/PixelDetector.cpp b/src/core/detector/PixelDetector.cpp index bc429e9c..7ca1bdc3 100644 --- a/src/core/detector/PixelDetector.cpp +++ b/src/core/detector/PixelDetector.cpp @@ -208,7 +208,7 @@ void PixelDetector::configure_pos_and_orientation(Configuration& config) const { PositionVector3D> PixelDetector::getIntercept(const Track* track) const { // FIXME: this is else statement can only be temporary - if(track->getType() == "GblTrack") { + if(track->getType() == "GblTrack" || track->getType() == "HelixTrack") { return track->getState(getName()); } else { // Get the distance from the plane to the track initial state diff --git a/src/objects/CMakeLists.txt b/src/objects/CMakeLists.txt index 246b7b6d..253a9902 100644 --- a/src/objects/CMakeLists.txt +++ b/src/objects/CMakeLists.txt @@ -27,6 +27,7 @@ ROOT_GENERATE_DICTIONARY(CorryvreckanObjectsDictionary ${CMAKE_CURRENT_SOURCE_DIR}/MCParticle.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Event.hpp ${CMAKE_CURRENT_SOURCE_DIR}/Multiplet.hpp + ${CMAKE_CURRENT_SOURCE_DIR}/HelixTrack.hpp LINKDEF ${CMAKE_CURRENT_SOURCE_DIR}/Linkdef.h OPTIONS @@ -52,6 +53,7 @@ ADD_LIBRARY(CorryvreckanObjects SHARED GblTrack.cpp Event.cpp Multiplet.cpp + HelixTrack.cpp ${CMAKE_CURRENT_BINARY_DIR}/CorryvreckanObjectsDictionary.cxx ) diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp new file mode 100644 index 00000000..7eab7384 --- /dev/null +++ b/src/objects/HelixTrack.cpp @@ -0,0 +1,383 @@ +/** + * @file + * @brief Implementation of StraightLine track object + * + * @copyright Copyright (c) 2017-2022 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#include "Eigen/Dense" +#include "StraightLineTrack.hpp" +#include "Track.hpp" +#include "core/utils/log.h" +#include "core/utils/unit.h" +#include "exceptions.h" + +using namespace corryvreckan; +using namespace Eigen; +using namespace ROOT::Math; + +ROOT::Math::XYPoint HelixTrack::getKinkAt(const std::string&) const { + return ROOT::Math::XYPoint(0, 0); +} + +ROOT::Math::XYZPoint HelixTrack::getState(const std::string& detectorID) const { + const Plane* plane = get_plane(detectorID); + XYZVector planeU, planeV, planeN; + plane->getToGlobal().Rotation().GetComponents(planeU, planeV, planeN); + return getIntercept(plane->getPosition(), planeN); +} + +ROOT::Math::XYZVector HelixTrack::getDirection(const std::string& detectorID) const { + const Plane* plane = get_plane(detectorID); + XYZVector planeU, planeV, planeN; + plane->getToGlobal().Rotation().GetComponents(planeU, planeV, planeN); + return getDirection(plane->getPosition(), planeN); +} + +ROOT::Math::XYZVector HelixTrack::getDirection(const double& z) const { + XYZPoint origin = XYZPoint(0, 0, z); + XYZVector normal = XYZVector(0, 0, 1); + return getDirection(origin, normal); +} + +void HelixTrack::calculate_transformations() { + // compute transformation to helix coordinate system with the direction of the magnetic field as the z-direction + // If there is no magnetic field, the transformation will simply result in identity + XYZVector globalZ = XYZVector(0, 0, 1); + XYZVector helixZ = magneticField_.R() != 0 ? magneticField_.unit() : globalZ; + XYZVector helixY = helixZ.Theta() != 0 ? helixZ.Cross(globalZ).unit() : XYZVector(0, 1, 0); + XYZVector helixX = helixY.Cross(helixZ).unit(); + + toGlobal_ = Rotation3D(helixX, helixY, helixZ); + toHelix_ = toGlobal_.Inverse(); +} + +void HelixTrack::calculate_chi2() { + + // Straightline: Each hit provides two measurements and there are four fitting parameters. + if(magneticField_.R() == 0) { + ndof_ = (track_clusters_.size() - 2) * 2; + } + // Helix: Each hit provides two measurements and there are five fitting paramters. + else { + ndof_ = 2 * track_clusters_.size() - 5; + } + chi2_ = 0.; + chi2ndof_ = 0.; + + // Loop over all clusters + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + // Get the distance and the error + auto intercept = getState(cluster->detectorID()); + auto interceptLocal = get_plane(cluster->detectorID())->getToLocal() * intercept; + auto residual = cluster->local() - interceptLocal; + + double ex2 = cluster->errorX() * cluster->errorX(); + double ey2 = cluster->errorY() * cluster->errorY(); + chi2_ += ((residual.x() * residual.x() / ex2) + (residual.y() * residual.y() / ey2)); + } + + // Store also the chi2/degrees of freedom + chi2ndof_ = (ndof_ <= 0) ? -1 : (chi2_ / static_cast(ndof_)); +} + +void HelixTrack::calculate_residuals() { + for(const auto& c : track_clusters_) { + auto* cluster = c.get(); + + auto intercept = getState(cluster->detectorID()); + auto interceptLocal = get_plane(cluster->detectorID())->getToLocal() * intercept; + residual_global_[cluster->detectorID()] = cluster->global() - intercept; + residual_local_[cluster->detectorID()] = cluster->local() - interceptLocal; + } +} + +double HelixTrack::operator()(const double* parameters) { + + // Update the HelixTrack parameters + this->curvature_ = parameters[0]; + this->phi_ = parameters[1]; + this->dca_ = parameters[2]; + this->tanLambda_ = parameters[3]; + this->z0_ = parameters[4]; + + // Calculate the chi2 + this->calculate_chi2(); + + // Return this to minuit + return chi2_; +} + +void HelixTrack::fit_helix() { + LOG(TRACE) << "Starting helix fit."; + // Fit points to helix with fast helix fit described in: + // https://www.physi.uni-heidelberg.de/Publications/DiplomaKiehn.pdf + + // 1. Fit circle in transverse plane + // + // mean of positons and radius in transverse plane (x, y, r^2) + Matrix mean = Matrix::Zero(); + // mean of squared positions and radius in transverse plane + // (x^2, x*y, x*r^2 ) + // ( 0 , y^2, y*r^2 ) + // ( 0 , 0 , r^4 ) + Matrix mean2 = Matrix::Zero(); + auto n = track_clusters_.size(); + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + auto position = toHelix_ * cluster->global(); + LOG(TRACE) << position; + long double x = position.x(); + long double y = position.y(); + long double r2 = x * x + y * y; + + mean(0) += x; + mean(1) += y; + mean(2) += r2; + mean2(0, 0) += x * x; + mean2(0, 1) += x * y; + mean2(0, 2) += x * r2; + mean2(1, 1) += y * y; + mean2(1, 2) += y * r2; + mean2(2, 2) += r2 * r2; + } + mean /= static_cast(n); + mean2 /= static_cast(n); + LOG(TRACE) << "mean: " << mean; + + Matrix C = Matrix::Zero(); + C(0, 0) = mean2(0, 0) - mean(0) * mean(0); + C(0, 1) = mean2(0, 1) - mean(0) * mean(1); + C(0, 2) = mean2(0, 2) - mean(0) * mean(2); + C(1, 1) = mean2(1, 1) - mean(1) * mean(1); + C(1, 2) = mean2(1, 2) - mean(1) * mean(2); + C(2, 2) = mean2(2, 2) - mean(2) * mean(2); + + LOG(TRACE) << "C: " << C; + + long double q1 = (C(2, 2) * C(0, 1)) - (C(0, 2) * C(1, 2)); + long double q2 = (C(2, 2) * (C(0, 0) - C(1, 1))) - (C(0, 2) * C(0, 2)) + (C(1, 2) * C(1, 2)); + LOG(TRACE) << "q1: " << q1 << ", q2: " << q2; + + phi_ = static_cast(0.5 * atan(2. * q1 / q2)); + long double kappa = (sin(phi_) * C(0, 2) - cos(phi_) * C(1, 2)) / C(2, 2); + long double delta = -kappa * mean(2) + sin(phi_) * mean(0) - cos(phi_) * mean(1); + + curvature_ = static_cast(2. * kappa / sqrt(1. - 4. * delta * kappa)); + dca_ = static_cast(2. * delta / (1. + sqrt(1. - 4. * delta * kappa))); + + LOG(TRACE) << "Circle fit finished. Radius: " << Units::display(1 / curvature_, {"m", "km"}) + << ", DCA: " << Units::display(dca_, {"mm", "um"}) << ", Phi: " << Units::display(phi_, "deg"); + + // 2. Calculate arc length in transverse plane + // 3. Fit in longitudinal (a,z)-plane + Matrix2d mat(Matrix2d::Zero()); + Vector2d vec(Vector2d::Zero()); + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + auto position = toHelix_ * cluster->global(); + double a = getArcLength(position.x(), position.y()); + + double z = position.z(); + + vec += Vector2d(z, a * z); + Matrix2d err; + err << 1., a, a, a * a; + mat += err; + } + + // Check for singularities. + if(fabs(mat.determinant()) < std::numeric_limits::epsilon()) { + throw TrackFitError(typeid(this), "Martix inversion in straight line fit failed"); + } + + // Get the longitudinal parameters + Eigen::Vector2d res = mat.inverse() * vec; + + z0_ = res[0]; + tanLambda_ = res[1]; + LOG(TRACE) << "Longitudinal fit finished. Z0: " << Units::display(z0_, {"mm", "cm"}) << ", tan(lambda): " << tanLambda_; +} + +void HelixTrack::fit_straightLine() { + LOG(TRACE) << "No Magnetic Field. Starting straightline fit."; + Matrix4d mat(Matrix4d::Zero()); + Vector4d vec(Vector4d::Zero()); + + for(auto& cl : track_clusters_) { + auto* cluster = cl.get(); + if(cluster == nullptr) { + throw MissingReferenceException(typeid(*this), typeid(Cluster)); + } + + double x = cluster->global().x(); + double y = cluster->global().y(); + double z = cluster->global().z(); + double ex2 = cluster->errorX() * cluster->errorX(); + double ey2 = cluster->errorY() * cluster->errorY(); + + Vector2d meas(x, y); + Matrix2d V; + V << ex2, 0., 0., ey2; + Matrix C; + C << 1., z, 0., 0., 0., 0., 1., z; + + vec += C.transpose() * V.inverse() * meas; + mat += C.transpose() * V.inverse() * C; + } + + // Check for singularities. + if(fabs(mat.determinant()) < std::numeric_limits::epsilon()) { + throw TrackFitError(typeid(this), "Martix inversion in straight line fit failed"); + } + // Get the StraightLineTrack parameters + Eigen::Vector4d res = mat.inverse() * vec; + + double bx = res[0]; + double by = res[2]; + double mx = res[1]; + double my = res[3]; + + tanLambda_ = 1. / sqrt(mx * mx + my * my); + phi_ = acos(mx * tanLambda_); + dca_ = tanLambda_ * (bx * my - by * mx); + z0_ = -tanLambda_ * tanLambda_ * (bx * mx + by * my); + curvature_ = 0.; + LOG(TRACE) << "Straightline fit finished. DCA: " << Units::display(dca_, {"mm", "um"}) + << ", Phi: " << Units::display(phi_, "deg") << ", z0: " << Units::display(z0_, "mm") + << ", tan(lambda): " << tanLambda_; +} + +void HelixTrack::fit() { + + isFitted_ = false; + calculate_transformations(); + + if(magneticField_.R() == 0.) { + fit_straightLine(); + } else { + fit_helix(); + } + + // Calculate the chi2 + this->calculate_chi2(); + this->calculate_residuals(); + isFitted_ = true; +} + +ROOT::Math::XYZPoint HelixTrack::getIntercept(double z) const { + XYZPoint origin = XYZPoint(0., 0., z); + XYZVector normal = XYZVector(0., 0., 1); + return getIntercept(origin, normal); +} + +double HelixTrack::getArcLength(double x, double y) const { + // Straightline + if(curvature_ == 0.) { + return cos(phi_) * x + sin(phi_) * y; + } + // Helix + // Transform to cooridnate system with origin in center of track circle + Vector2d posOrigin; + posOrigin << x, y; + Matrix2d toCenter; + toCenter << sin(phi_), -cos(phi_), cos(phi_), sin(phi_); + Vector2d posCenter = toCenter * posOrigin; + posCenter[0] -= (dca_ + 1. / curvature_); + + double a = abs(1. / curvature_) * + (curvature_ <= 0. ? std::atan2(posCenter[1], posCenter[0]) : std::atan2(posCenter[1], -posCenter[0])); + return a; +} + +ROOT::Math::XYZPoint HelixTrack::getIntercept(XYZPoint origin, XYZVector normal) const { + // Rotate plane origin and normal to helix coordinates + auto normalHelix = toHelix_ * normal; + auto originHelix = toHelix_ * origin; + + double cosLambda = 1. / sqrt(1. + tanLambda_ * tanLambda_); + double sinLambda = sqrt(1. - cosLambda * cosLambda); + XYZVector distance, trackDirection; + XYZPoint position; + if(curvature_ == 0.) { + trackDirection = XYZVector(cosLambda * cos(phi_), cosLambda * sin(phi_), sinLambda); + // Point of closest approach + XYZPoint pca(dca_ * sin(phi_), -dca_ * cos(phi_), z0_); + distance = pca - originHelix; + double pathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + position = pca + pathLength * trackDirection; + } else { + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(dPhi - phi_); + const double sinPhi = sin(dPhi - phi_); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) + cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; + } + } + } + return toGlobal_ * position; +} + +ROOT::Math::XYZVector HelixTrack::getDirection(XYZPoint origin, XYZVector normal) const { + // Rotate plane normal vector to helix coordinates + auto normalHelix = toHelix_ * normal; + auto originHelix = toHelix_ * origin; + + double cosLambda = 1. / sqrt(1. + tanLambda_ * tanLambda_); + double sinLambda = sqrt(1. - cosLambda * cosLambda); + XYZVector distance, trackDirection; + XYZPoint position; + if(curvature_ == 0.) { + trackDirection = XYZVector(cosLambda * cos(phi_), cosLambda * sin(phi_), sinLambda); + } else { + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(phi_ + dPhi); + const double sinPhi = sin(phi_ + dPhi); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) - cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; + } + } + } + return toGlobal_ * trackDirection; +} + +void HelixTrack::print(std::ostream& out) const { + out << "HelixTrack " << this->curvature_ << ", " << this->phi_ << ", " << this->dca_ << ", " << this->tanLambda_ << ", " + << this->z0_ << ", " << this->chi2_ << ", " << this->ndof_ << ", " << this->chi2ndof_ << ", " << this->timestamp(); +} diff --git a/src/objects/HelixTrack.hpp b/src/objects/HelixTrack.hpp new file mode 100644 index 00000000..f7f84c5f --- /dev/null +++ b/src/objects/HelixTrack.hpp @@ -0,0 +1,125 @@ +/** + * @file + * @brief Definition of Helix track object + * + * @copyright Copyright (c) 2017-2020 CERN and the Corryvreckan authors. + * This software is distributed under the terms of the MIT License, copied verbatim in the file "LICENSE.md". + * In applying this license, CERN does not waive the privileges and immunities granted to it by virtue of its status as an + * Intergovernmental Organization or submit itself to any jurisdiction. + */ + +#ifndef CORRYVRECKAN_HELIXTRACK_H +#define CORRYVRECKAN_HELIXTRACK_H 1 + +#include "Track.hpp" + +namespace corryvreckan { + /** + * @ingroup Objects + * @brief HelixTrack object + * + * This class is a simple track class which knows how to fit itself. It is derived from Track + */ + + class HelixTrack : public Track { + + public: + // Minimisation operator used by Minuit. Minuit passes the current iteration of the parameters and checks if the chi2 + // is better or worse + double operator()(const double* parameters); + + void print(std::ostream& out) const override; + + /** + * @brief The fiting routine + */ + void fit() override; + + /** + * @brief Get the track position for a certain z position + * @param z Global z position + * @return ROOT::Math::XYZPoint at z position + */ + ROOT::Math::XYZPoint getIntercept(double z) const override; + + /** + * @brief Get the track position on a plane specified by a point and the normal vector in global coordinates + * @param origin Point on the plane (in global coordinates) + * @param normal Normal vector to the plane (in global coordinates) + * @return ROOT::Math::XYZPoint at intersection with plane + */ + ROOT::Math::XYZPoint getIntercept(ROOT::Math::XYZPoint origin, ROOT::Math::XYZVector normal) const; + + /** + * @brief Get the track state at a detector + * @param detectorID Name of detector + * @return ROOT::Math::XYZPoint state at detetcor layer + */ + ROOT::Math::XYZPoint getState(const std::string& detectorID) const override; + + /** + * @brief Get the track direction at a detector + * @param detectorID Name of detector + * @return ROOT::Math::XYZPoint direction at detetcor layer + */ + ROOT::Math::XYZVector getDirection(const std::string& detectorID) const override; + + /** + * @brief Get the track direction at any given z position + * @param z Global z position + * @return ROOT::Math::XYZPoint direction at z + */ + ROOT::Math::XYZVector getDirection(const double& z) const override; + + /** + * @brief Get the track direction on a plane specified by a point and the normal vector in global coordinates + * @param origin Point on the plane (in global coordinates) + * @param normal Normal vector to the plane (in global coordinates) + * @return ROOT::Math::XYZVector at intersection with plane + */ + ROOT::Math::XYZVector getDirection(ROOT::Math::XYZPoint origin, ROOT::Math::XYZVector normal) const; + + /** + * @brief This track model does not support kinks, it therefore is always zero + * @param detectorID Detector ID at which the kink should be evaluated + * @return Kink at given detector + */ + ROOT::Math::XYPoint getKinkAt(const std::string& detectorID) const override; + + /** + * @brief Get the arc length in the transverse plane of the helix + * @param x x-coordinate in transverse plane + * @param y y-coordinate in transverse plane + * @return double arc length + */ + double getArcLength(double x, double y) const; + + void setVolumeScatter(double) override{}; + + private: + /** + * @brief calculate the chi2 of the linear regression + */ + void calculate_chi2(); + + void calculate_residuals(); + + void calculate_transformations(); + + void fit_straightLine(); + void fit_helix(); + // Member variables (Karimaeki parametrization) + double curvature_; + double phi_; + double dca_; + double tanLambda_; + double z0_; + ROOT::Math::Rotation3D toHelix_; + ROOT::Math::Rotation3D toGlobal_; + + // ROOT I/O class definition - update version number when you change this class! + ClassDefOverride(HelixTrack, 1) + }; +} // namespace corryvreckan + +#endif // CORRYVRECKAN_HELIXTRACK_H diff --git a/src/objects/Linkdef.h b/src/objects/Linkdef.h index 9c35e804..04bf5730 100644 --- a/src/objects/Linkdef.h +++ b/src/objects/Linkdef.h @@ -23,6 +23,7 @@ #pragma link C++ class corryvreckan::StraightLineTrack + ; #pragma link C++ class corryvreckan::GblTrack + ; #pragma link C++ class corryvreckan::Multiplet + ; +#pragma link C++ class corryvreckan::HelixTrack + ; #pragma link C++ class corryvreckan::MCParticle + ; #pragma link C++ class corryvreckan::Event + ; #pragma link C++ class corryvreckan::Track::Plane + ; diff --git a/src/objects/Track.cpp b/src/objects/Track.cpp index 11932f1d..1c238f00 100644 --- a/src/objects/Track.cpp +++ b/src/objects/Track.cpp @@ -240,7 +240,7 @@ void Track::registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, } } -Track::Plane* Track::get_plane(std::string detetorID) { +const Track::Plane* Track::get_plane(const std::string detetorID) const { auto plane = std::find_if(planes_.begin(), planes_.end(), [&detetorID](Plane const& p) { return p.getName() == detetorID; }); if(plane == planes_.end()) { @@ -254,6 +254,8 @@ std::shared_ptr corryvreckan::Track::Factory(const std::string& trackMode return std::make_shared(); } else if(trackModel == "gbl") { return std::make_shared(); + } else if(trackModel == "helix") { + return std::make_shared(); } else { throw UnknownTrackModel(typeid(Track), trackModel); } diff --git a/src/objects/Track.hpp b/src/objects/Track.hpp index 9b799ed5..1a2a5148 100644 --- a/src/objects/Track.hpp +++ b/src/objects/Track.hpp @@ -299,7 +299,7 @@ namespace corryvreckan { void loadHistory() override; void petrifyHistory() override; - Plane* get_plane(std::string detetorID); + const Plane* get_plane(const std::string detetorID) const; std::vector> track_clusters_; std::vector> associated_clusters_; std::map residual_local_; @@ -324,6 +324,7 @@ namespace corryvreckan { // include all tracking methods here to have one header to be include everywhere #include "GblTrack.hpp" +#include "HelixTrack.hpp" #include "Multiplet.hpp" #include "StraightLineTrack.hpp" #endif // CORRYVRECKAN_TRACK_H diff --git a/src/objects/objects.h b/src/objects/objects.h index ff247b84..051336e4 100644 --- a/src/objects/objects.h +++ b/src/objects/objects.h @@ -18,5 +18,5 @@ namespace corryvreckan { /** * @brief Tuple containing all objects */ - using OBJECTS = std::tuple; + using OBJECTS = std::tuple; } // namespace corryvreckan -- GitLab From 53c1184c50d3c9ffd649596b523ac800d52dc60d Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 21 Jun 2022 11:56:12 +0200 Subject: [PATCH 05/14] Tracking4D: Use helix as reference track --- src/modules/Tracking4D/Tracking4D.cpp | 21 +++++++++++++++++---- 1 file changed, 17 insertions(+), 4 deletions(-) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index d764ed08..b37cdd24 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -24,10 +24,11 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("min_hits_on_track", 6); config_.setDefault("exclude_dut", true); - config_.setDefault("track_model", "straightline"); + config_.setDefault("track_model", "helix"); config_.setDefault("momentum", Units::get(5, "GeV")); config_.setDefault("max_plot_chi2", 50.0); config_.setDefault("volume_radiation_length", Units::get(304.2, "m")); @@ -308,9 +309,17 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { } // The track finding is based on a straight line. Therefore a refTrack to extrapolate to the next plane is used - StraightLineTrack refTrack; + HelixTrack refTrack; refTrack.addCluster(clusterFirst.get()); + refTrack.registerPlane(reference_first->getName(), + reference_first->origin(), + reference_first->materialBudget(), + reference_first->toLocal()); refTrack.addCluster(clusterLast.get()); + refTrack.registerPlane(reference_last->getName(), + reference_last->origin(), + reference_last->materialBudget(), + reference_last->toLocal()); auto averageTimestamp = calculate_average_timestamp(&refTrack); refTrack.setTimestamp(averageTimestamp); @@ -323,6 +332,8 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { if(use_volume_scatterer_) { track->setVolumeScatter(volume_radiation_length_); } + refTrack.setParticleMomentum(momentum_); + refTrack.setMagneticField(XYZVector(0, 0, 0)); track->setParticleMomentum(momentum_); track->setMagneticField(magnetic_field_); @@ -335,8 +346,8 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { } auto detectorID = detector->getName(); LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); - track->registerPlane( - detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); + refTrack.registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); + track->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); if(detector == reference_first || detector == reference_last) { continue; @@ -432,6 +443,8 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { averageTimestamp = calculate_average_timestamp(&refTrack); refTrack.setTimestamp(averageTimestamp); track->setTimestamp(averageTimestamp); + // At least 3 clusters, helix fit is possible + refTrack.setMagneticField(magnetic_field_); LOG(DEBUG) << "- added cluster to track"; } -- GitLab From d1582e066bd0b05ed951f19046a1c9c633fca9aa Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 21 Jun 2022 13:43:02 +0200 Subject: [PATCH 06/14] Add curvature from magnetic field to GblTrack --- src/objects/GblTrack.cpp | 37 ++++++++++++++++++++++++------------- 1 file changed, 24 insertions(+), 13 deletions(-) diff --git a/src/objects/GblTrack.cpp b/src/objects/GblTrack.cpp index a37c84e1..518c72c8 100644 --- a/src/objects/GblTrack.cpp +++ b/src/objects/GblTrack.cpp @@ -56,6 +56,10 @@ void GblTrack::add_plane(std::vector::iterator& plane, double total_material) { // lambda to add plane (not the first one) and air scatterers auto globalTangent = Vector4d(0, 0, 1, 0); + auto cosLamb = std::sqrt(globalTangent.x() * globalTangent.x() + globalTangent.y() * globalTangent.y()) / + std::sqrt(globalTangent.x() * globalTangent.x() + globalTangent.y() * globalTangent.y() + + globalTangent.z() * globalTangent.z()); + auto lamb = std::acos(cosLamb); // extract the rotation from an ROOT::Math::Transfrom3D, store it in 4x4 matrix to match proteus format auto getRotation = [](Transform3D in) { @@ -69,13 +73,13 @@ void GblTrack::add_plane(std::vector::iterator& plane, auto tmp_local = plane->getToLocal() * globalTrackPos; Vector4d localPosTrack = Vector4d{tmp_local.x(), tmp_local.y(), tmp_local.z(), 1}; Vector4d localTangent = getRotation(plane->getToLocal()) * globalTangent; - double dist = localPosTrack[2]; + double dist = -localPosTrack[2]; // Minus sign to get a positive value LOG(TRACE) << "Rotation: " << getRotation(plane->getToLocal()); LOG(TRACE) << "Local tan before normalization: " << localTangent; LOG(TRACE) << "Distance: " << dist; localTangent /= localTangent.z(); - localPosTrack -= dist * localTangent; + localPosTrack += dist * localTangent; // add the local track pos for future reference - e.g. dut position: local_track_points_[plane->getName()] = ROOT::Math::XYPoint(localPosTrack(0), localPosTrack(1)); @@ -83,11 +87,15 @@ void GblTrack::add_plane(std::vector::iterator& plane, Matrix4d toTarget = getRotation(plane->getToLocal()) * getRotation(prevToGlobal); // lambda Jacobian from one scatter to the next - auto jac = [](const Vector4d& tangent, const Matrix4d& target, double distance) { + auto jac = [this](const Vector4d& tangent, const Matrix4d& target, double distance, double lambda) { + auto magneticFieldFactor = magneticField_ * 0.3; // Magnetic field factor: B*c/MeV*mm + auto BfacX = (-std::cos(lambda) * magneticFieldFactor.z() - std::sin(lambda) * magneticFieldFactor.y()) * distance; + auto BfacY = magneticFieldFactor.x() * distance; + LOG(TRACE) << "Magnetic field factor: " << magneticFieldFactor; Matrix R; R.col(0) = target.col(0); R.col(1) = target.col(1); - R.col(2) = target.col(3); + R.col(2) = target.col(0) * BfacX + target.col(1) * BfacY; Vector4d S = target * tangent * (1 / tangent[2]); Matrix F = Matrix::Zero(); @@ -99,11 +107,13 @@ void GblTrack::add_plane(std::vector::iterator& plane, F(2, 2) = -S[3] / S[2]; Matrix jaco; - jaco << F * R, (-distance / S[2]) * F * R, Matrix3d::Zero(), (1 / S[2]) * F * R; - jaco(5, 5) = 1; // a future time component + jaco << F * R, (distance / S[2]) * F * R, Matrix3d::Zero(), (1 / S[2]) * F * R; + jaco(0, 5) /= 2; + jaco(1, 5) /= 2; + jaco(5, 5) = 1; return jaco; }; - auto myjac = jac(prevTan, toTarget, dist); + auto myjac = jac(prevTan, toTarget, dist, lamb); // Layout if volume scattering active // | | | | @@ -142,16 +152,16 @@ void GblTrack::add_plane(std::vector::iterator& plane, myjac(0, 0) = 1; // Adding volume scattering if requested } else if(use_volume_scatter_) { - myjac = jac(prevTan, toTarget, frac1 * dist); + myjac = jac(prevTan, toTarget, frac1 * dist, lamb); GblPoint pVolume(toGbl * myjac * toProt); addScattertoGblPoint(pVolume, fabs(dist) / 2. / scattering_length_volume_); gblpoints_.push_back(pVolume); // We have already rotated to the next local coordinate system - myjac = jac(prevTan, Matrix4d::Identity(), frac2 * dist); + myjac = jac(prevTan, Matrix4d::Identity(), frac2 * dist, lamb); GblPoint pVolume2(toGbl * myjac * toProt); addScattertoGblPoint(pVolume2, fabs(dist) / 2. / scattering_length_volume_); gblpoints_.push_back(pVolume2); - myjac = jac(prevTan, Matrix4d::Identity(), frac1 * dist); + myjac = jac(prevTan, Matrix4d::Identity(), frac1 * dist, lamb); } auto transformedJac = toGbl * myjac * toProt; GblPoint point(transformedJac); @@ -267,7 +277,7 @@ void GblTrack::fit() { prepare_gblpoints(); // perform fit - GblTrajectory traj(gblpoints_, false); // false = no magnetic field + GblTrajectory traj(gblpoints_, magneticField_.R() != 0); // true = magnetic field, false = no magnetic field double lostWeight = 0; int ndf = 0; @@ -337,7 +347,7 @@ void GblTrack::fit() { LOG(TRACE) << "Results for detector " << name << std::endl << "Fitted residual local:\t" << residual_local_.at(name) << std::endl << "Seed residual:\t" << initital_residual_.at(name) << std::endl - << "Ditted residual global:\t" << ROOT::Math::XYPoint(clusterPos - corPos); + << "Fitted residual global:\t" << ROOT::Math::XYPoint(clusterPos - corPos); } LOG(DEBUG) << "Plane: " << name << ": residual " << residual_local_[name] << ", kink: " << kink_[name]; } @@ -447,7 +457,8 @@ ROOT::Math::XYZVector GblTrack::getDirection(const std::string& detectorID) cons } XYZVector GblTrack::getDirection(const double& z) const { - auto planeUpstream = std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition().z() > z; }); + auto planeUpstream = + std::find_if(planes_.begin(), planes_.end(), [&z](const Plane& p) { return p.getPosition().z() > z; }); if(planeUpstream != planes_.end()) { return getDirection(planeUpstream->getName()); } else { -- GitLab From 9aa7639623241c4b552f4f813615aed9b26569f4 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 21 Jun 2022 17:48:47 +0200 Subject: [PATCH 07/14] HelixTrack: Handle straightline case with zero slope --- src/modules/Tracking4D/Tracking4D.cpp | 2 +- src/objects/HelixTrack.cpp | 20 ++++++++++++++------ 2 files changed, 15 insertions(+), 7 deletions(-) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index b37cdd24..aaf2b559 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -24,7 +24,7 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("min_hits_on_track", 6); config_.setDefault("exclude_dut", true); diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp index 7eab7384..1e1779f6 100644 --- a/src/objects/HelixTrack.cpp +++ b/src/objects/HelixTrack.cpp @@ -254,13 +254,21 @@ void HelixTrack::fit_straightLine() { double mx = res[1]; double my = res[3]; - tanLambda_ = 1. / sqrt(mx * mx + my * my); - phi_ = acos(mx * tanLambda_); - dca_ = tanLambda_ * (bx * my - by * mx); - z0_ = -tanLambda_ * tanLambda_ * (bx * mx + by * my); + // cover case without any slope in x or y direction + if(mx < std::numeric_limits::epsilon() && my < std::numeric_limits::epsilon()) { + tanLambda_ = std::numeric_limits::max(); + phi_ = atan2(bx, -by); + dca_ = sin(phi_) * bx - cos(phi_) * by; + z0_ = 0.; + } else { + tanLambda_ = 1. / sqrt(mx * mx + my * my); + phi_ = acos(mx * tanLambda_); + dca_ = tanLambda_ * (bx * my - by * mx); + z0_ = -(bx * mx + by * my) * tanLambda_ * tanLambda_; + } curvature_ = 0.; - LOG(TRACE) << "Straightline fit finished. DCA: " << Units::display(dca_, {"mm", "um"}) - << ", Phi: " << Units::display(phi_, "deg") << ", z0: " << Units::display(z0_, "mm") + LOG(TRACE) << "Straightline fit finished. DCA: " << Units::display(dca_, {"m", "cm", "mm", "um"}) + << ", Phi: " << Units::display(phi_, "deg") << ", z0: " << Units::display(z0_, {"m", "cm", "mm", "um"}) << ", tan(lambda): " << tanLambda_; } -- GitLab From 6916b56589185fc8e5cbd48157aa7e0a35821343 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Mon, 11 Jul 2022 15:13:57 +0200 Subject: [PATCH 08/14] Fix calculation of phi for straight line case --- src/objects/HelixTrack.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp index 1e1779f6..58d8964b 100644 --- a/src/objects/HelixTrack.cpp +++ b/src/objects/HelixTrack.cpp @@ -79,6 +79,7 @@ void HelixTrack::calculate_chi2() { auto intercept = getState(cluster->detectorID()); auto interceptLocal = get_plane(cluster->detectorID())->getToLocal() * intercept; auto residual = cluster->local() - interceptLocal; + LOG(TRACE) << "Residual on detector " << cluster->detectorID() << ": " << residual * 1000 << " um"; double ex2 = cluster->errorX() * cluster->errorX(); double ey2 = cluster->errorY() * cluster->errorY(); @@ -254,6 +255,8 @@ void HelixTrack::fit_straightLine() { double mx = res[1]; double my = res[3]; + LOG(TRACE) << "bx = " << bx << ", by = " << by << ", mx = " << mx << ", my = " << my; + // cover case without any slope in x or y direction if(mx < std::numeric_limits::epsilon() && my < std::numeric_limits::epsilon()) { tanLambda_ = std::numeric_limits::max(); @@ -262,7 +265,7 @@ void HelixTrack::fit_straightLine() { z0_ = 0.; } else { tanLambda_ = 1. / sqrt(mx * mx + my * my); - phi_ = acos(mx * tanLambda_); + phi_ = atan2(my, mx); dca_ = tanLambda_ * (bx * my - by * mx); z0_ = -(bx * mx + by * my) * tanLambda_ * tanLambda_; } -- GitLab From 7416d8395c5c2b5e5b84081858730ffae6fe8bae Mon Sep 17 00:00:00 2001 From: David Bacher Date: Mon, 11 Jul 2022 15:15:50 +0200 Subject: [PATCH 09/14] Add magnetic field to Prealignment --- src/modules/Prealignment/Prealignment.cpp | 47 +++++++++++++++++++++++ src/modules/Prealignment/Prealignment.h | 5 +++ 2 files changed, 52 insertions(+) diff --git a/src/modules/Prealignment/Prealignment.cpp b/src/modules/Prealignment/Prealignment.cpp index 262ba4b9..de8cdc64 100644 --- a/src/modules/Prealignment/Prealignment.cpp +++ b/src/modules/Prealignment/Prealignment.cpp @@ -9,6 +9,7 @@ */ #include "Prealignment.h" +#include #include "tools/cuts.h" using namespace corryvreckan; @@ -25,6 +26,10 @@ Prealignment::Prealignment(Configuration& config, std::shared_ptr dete config_.setDefault("method", PrealignMethod::MEAN); config_.setDefault("fit_range_rel", 500); config_.setDefault("range_abs", Units::get(10, "mm")); + config_.setDefault("magnetic_field", ROOT::Math::XYZVector()); + config_.setDefault("momentum", Units::get(5, "GeV")); + config_.setDefault("charge", Units::get(1, "e")); + config_.setDefault("initial_angle", Units::get(0, "deg")); if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { config_.setDefault("time_cut_rel", 3.0); @@ -38,6 +43,10 @@ Prealignment::Prealignment(Configuration& config, std::shared_ptr dete range_abs = config_.get("range_abs"); method = config_.get("method"); fit_range_rel = config_.get("fit_range_rel"); + magnetic_field = config_.get("magnetic_field"); + momentum = config_.get("momentum"); + charge = config_.get("charge"); + initial_angle = config_.get("initial_angle"); LOG(DEBUG) << "Setting max_correlation_rms to : " << max_correlation_rms; LOG(DEBUG) << "Setting damping_factor to : " << damping_factor; @@ -167,6 +176,44 @@ void Prealignment::finalize(const std::shared_ptr&) { shift_Y = correlationY->GetXaxis()->GetBinCenter(binMaxY); } + if(magnetic_field.R() != 0) { + ROOT::Math::XYZVector globalZ(0, 0, 1); + auto helixZ = magnetic_field.unit(); + auto helixY = helixZ.Theta() != 0 ? helixZ.Cross(globalZ).unit() : ROOT::Math::XYZVector(0, 1, 0); + auto helixX = helixY.Cross(helixZ).unit(); + + ROOT::Math::Rotation3D toGlobal(helixX, helixY, helixZ); + auto toHelix = toGlobal.Inverse(); + + // calculate bending radius (in mm) from magnetic field, momentum and charge + double R = momentum / charge / magnetic_field.R() * globalZ.Cross(magnetic_field.unit()).R() * 10. / 3.; + + ROOT::Math::XYZVector planePosGlobal(0., 0., m_detector->displacement().Z()); + ROOT::Math::XYZVector refPosGlobal(0., 0., get_reference()->displacement().Z()); + auto planePosHelix = toHelix * planePosGlobal; + auto refPosHelix = toHelix * refPosGlobal; + planePosHelix.SetX(planePosHelix.x() + R * sin(initial_angle)); + refPosHelix.SetX(refPosHelix.x() + R * sin(initial_angle)); + + double initialShift = -(charge / abs(charge)) * R * (1 - cos(initial_angle)); + double planeShiftHelix = + -(charge / abs(charge)) * (R - sqrt(R * R - planePosHelix.x() * planePosHelix.x())) - initialShift; + double refShiftHelix = + -(charge / abs(charge)) * (R - sqrt(R * R - refPosHelix.x() * refPosHelix.x())) - initialShift; + auto planeShiftGlobal = toGlobal * ROOT::Math::XYZVector(0., planeShiftHelix, 0.); + auto refShiftGlobal = toGlobal * ROOT::Math::XYZVector(0., refShiftHelix, 0.); + + auto magneticShift_X = planeShiftGlobal.x() - refShiftGlobal.x(); + auto magneticShift_Y = planeShiftGlobal.y() - refShiftGlobal.y(); + + shift_X += magneticShift_X; + shift_Y += magneticShift_Y; + + LOG(INFO) << "Shift (from magnetic field)" << m_detector->getName() + << ": x = " << Units::display(magneticShift_X, {"mm", "um"}) + << " , y = " << Units::display(magneticShift_Y, {"mm", "um"}); + } + LOG(DEBUG) << "Shift (without damping factor)" << m_detector->getName() << ": x = " << Units::display(shift_X, {"mm", "um"}) << " , y = " << Units::display(shift_Y, {"mm", "um"}); diff --git a/src/modules/Prealignment/Prealignment.h b/src/modules/Prealignment/Prealignment.h index 4a584f97..9ebc935a 100644 --- a/src/modules/Prealignment/Prealignment.h +++ b/src/modules/Prealignment/Prealignment.h @@ -11,6 +11,7 @@ #ifndef PREALIGNMENT_H #define PREALIGNMENT_H 1 +#include #include #include #include @@ -59,6 +60,10 @@ namespace corryvreckan { double range_abs; PrealignMethod method; int fit_range_rel; + ROOT::Math::XYZVector magnetic_field; + double momentum; + double charge; + double initial_angle; }; } // namespace corryvreckan #endif // PREALIGNMENT_H -- GitLab From 2578ad6d56c738d85c1dd579342d9f73e4b42fad Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 12 Jul 2022 15:29:23 +0200 Subject: [PATCH 10/14] Throw error if straightlines are used with magnetic field --- src/modules/Tracking4D/Tracking4D.cpp | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index aaf2b559..b89532ee 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -70,6 +70,12 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("unique_cluster_usage"); magnetic_field_ = config_.get("magnetic_field"); + // straightline tracks with magnetic field are unphysical + if(track_model_ == "straightline" && magnetic_field_ != ROOT::Math::XYZVector()) { + throw InvalidCombinationError( + config_, {"track_model", "magnetic_field"}, "Straightlines do not support magnetic fields."); + } + // print a warning if volumeScatterer are used as this causes fit failures // that are still not understood if(use_volume_scatterer_) { -- GitLab From d68b2e5b4d4d6400f112ef7ec8393d449b2f95d9 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Wed, 13 Jul 2022 17:50:17 +0200 Subject: [PATCH 11/14] Tracking4D: track seeding with different models `StraightLine` if no magnetic field, `Helix` with magnetic field --- src/modules/Tracking4D/Tracking4D.cpp | 405 +++++++++++++++----------- 1 file changed, 231 insertions(+), 174 deletions(-) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index cb36fcc1..0f4e82db 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -271,7 +271,16 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { // Container for all clusters, and detectors in tracking map, KDTree> trees; - std::shared_ptr reference_first, reference_last; + std::vector> seed_detectors; + // Number of detectors needed for seeding (2 without magnetic field, 3 with magnetic field) + size_t num_seed_detectors; + if(magnetic_field_ == ROOT::Math::XYZVector()) { + LOG(DEBUG) << "No magnetic field. Seeding tracks from two detectors."; + num_seed_detectors = 2; + } else { + LOG(DEBUG) << "Mangetic field. Seeding track from three detectors."; + num_seed_detectors = 3; + } for(auto& detector : get_regular_detectors(!exclude_DUT_)) { // Get the clusters auto tempClusters = clipboard->getData(detector->getName()); @@ -286,10 +295,11 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { // Get first and last detectors with clusters on them: if(std::find(exclude_from_seed_.begin(), exclude_from_seed_.end(), detector->getName()) == exclude_from_seed_.end()) { - if(!reference_first) { - reference_first = detector; + if(seed_detectors.size() < num_seed_detectors) { + seed_detectors.push_back(detector); + } else { + seed_detectors.back() = detector; } - reference_last = detector; } else { LOG(DEBUG) << "Not using " << detector->getName() << " as seed as chosen by config file."; } @@ -297,7 +307,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { } // If there are no detectors then stop trying to track - if(trees.size() < 2) { + if(seed_detectors.size() < num_seed_detectors) { // Fill histogram tracksPerEvent->Fill(0); @@ -309,213 +319,260 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { TrackVector tracks; // Time cut for combinations of reference clusters and for reference track with additional detector - auto time_cut_ref = std::max(time_cuts_[reference_first], time_cuts_[reference_last]); - auto time_cut_ref_track = std::min(time_cuts_[reference_first], time_cuts_[reference_last]); - for(auto& clusterFirst : trees[reference_first].getAllElements()) { - for(auto& clusterLast : trees[reference_last].getAllElements()) { - LOG(DEBUG) << "Looking at next reference cluster pair"; - - if(std::fabs(clusterFirst->timestamp() - clusterLast->timestamp()) > time_cut_ref) { - LOG(DEBUG) << "Reference clusters not within time cuts."; - continue; + double time_cut_ref = std::numeric_limits::min(); + double time_cut_ref_track = std::numeric_limits::max(); + std::map, std::vector>> seed_clusters; + std::map, std::vector>::iterator> seed_clusters_it; + for(auto seedDetector : seed_detectors) { + time_cut_ref = std::max(time_cut_ref, time_cuts_[seedDetector]); + time_cut_ref_track = std::min(time_cut_ref_track, time_cuts_[seedDetector]); + seed_clusters[seedDetector] = trees[seedDetector].getAllElements(); + seed_clusters_it[seedDetector] = seed_clusters[seedDetector].begin(); + } + + // lambda to iterate through all combinations of elements of variable number of vectors + auto iterate = [&seed_detectors, &seed_clusters, &seed_clusters_it]() { + ++seed_clusters_it[seed_detectors.front()]; + for(auto seedDetectorsIt = seed_detectors.begin(); seedDetectorsIt != seed_detectors.end() - 1; seedDetectorsIt++) { + if(seed_clusters_it[*seedDetectorsIt] != seed_clusters[*seedDetectorsIt].end()) { + break; } + seed_clusters_it[*seedDetectorsIt] = seed_clusters[*seedDetectorsIt].begin(); + ++seed_clusters_it[*(seedDetectorsIt + 1)]; + } + }; - // The track finding is based on a straight line. Therefore a refTrack to extrapolate to the next plane is used - HelixTrack refTrack; - refTrack.addCluster(clusterFirst.get()); - refTrack.registerPlane(reference_first->getName(), - reference_first->origin(), - reference_first->materialBudget(), - reference_first->toLocal()); - refTrack.addCluster(clusterLast.get()); - refTrack.registerPlane(reference_last->getName(), - reference_last->origin(), - reference_last->materialBudget(), - reference_last->toLocal()); - auto averageTimestamp = calculate_average_timestamp(&refTrack); - refTrack.setTimestamp(averageTimestamp); - - // Make a new track - auto track = Track::Factory(track_model_); - track->addCluster(clusterFirst.get()); - track->addCluster(clusterLast.get()); + for(; seed_clusters_it[seed_detectors.back()] != seed_clusters[seed_detectors.back()].end(); iterate()) { + LOG(DEBUG) << "Looking at next reference cluster combination"; - track->setTimestamp(averageTimestamp); - if(use_volume_scatterer_) { - track->setVolumeScatter(volume_radiation_length_); - } - refTrack.setParticleMomentum(momentum_); - refTrack.setMagneticField(XYZVector(0, 0, 0)); - track->setParticleMomentum(momentum_); - track->setMagneticField(magnetic_field_); - - // Loop over each subsequent plane and look for a cluster within the timing cuts - size_t detector_nr = 2; - // Get all detectors here to also include passive layers which might contribute to scattering - for(auto& detector : get_detectors()) { - if(detector->isAuxiliary()) { - continue; - } - auto detectorID = detector->getName(); - LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); - refTrack.registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); - track->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); + std::map, std::shared_ptr> seedClusters; + for(auto& seedDetector : seed_detectors) { + seedClusters[seedDetector] = *seed_clusters_it[seedDetector]; + } - if(detector == reference_first || detector == reference_last) { - continue; + auto checkSeedTimeCuts = [&seed_detectors, + time_cut_ref](std::map, std::shared_ptr>& clusters) { + auto seedDetector1 = seed_detectors.begin(); + while(seedDetector1 != seed_detectors.end()) { + auto seedDetector2 = seedDetector1 + 1; + while(seedDetector2 != seed_detectors.end()) { + if(std::fabs(clusters[*seedDetector1]->timestamp() - clusters[*seedDetector2]->timestamp()) > + time_cut_ref) { + LOG(DEBUG) << "Reference clusters on " << (*seedDetector1)->getName() << " and " + << (*seedDetector2)->getName() << " not within time cuts."; + return false; + } + seedDetector2++; } + seedDetector1++; + } + return true; + }; - if(exclude_DUT_ && detector->isDUT()) { - LOG(DEBUG) << "Skipping DUT plane."; - continue; - } + if(!checkSeedTimeCuts(seedClusters)) { + continue; + } - if(detector->isPassive()) { - LOG(DEBUG) << "Skipping passive plane."; - continue; - } + // The track finding is based on a straight line or a helix. Therefore a seedTrack to extrapolate to the next plane + // is used + std::string seedTrackModel; + if(magnetic_field_ == ROOT::Math::XYZVector()) { + LOG(DEBUG) << "Seeding tracks based on straight line model."; + seedTrackModel = "straightline"; + } else { + LOG(DEBUG) << "Seeding tracks based on helix model."; + seedTrackModel = "helix"; + } - // Determine whether a track can still be assembled given the number of current hits and the number of - // detectors to come. Reduces computing time. - detector_nr++; - if(refTrack.getNClusters() + (trees.size() - detector_nr + 1) < min_hits_on_track_) { - LOG(DEBUG) << "No chance to find a track - too few detectors left: " << refTrack.getNClusters() << " + " - << trees.size() << " - " << detector_nr << " < " << min_hits_on_track_; - continue; - } + // Make new seed track and new track + auto seedTrack = Track::Factory(seedTrackModel); + auto track = Track::Factory(track_model_); - if(trees.count(detector) == 0) { - LOG(TRACE) << "Skipping detector " << detector->getName() << " as it has 0 clusters."; - continue; - } + for(auto& seedDetector : seed_detectors) { + seedTrack->addCluster(seedClusters[seedDetector].get()); + seedTrack->registerPlane( + seedDetector->getName(), seedDetector->origin(), seedDetector->materialBudget(), seedDetector->toLocal()); - // Get all neighbors within the timing cut - LOG(DEBUG) << "Searching for neighboring cluster on device " << detector->getName(); - LOG(DEBUG) << "- reference time is " << Units::display(refTrack.timestamp(), {"ns", "us", "s"}); - Cluster* closestCluster = nullptr; + track->addCluster(seedClusters[seedDetector].get()); + } - // Use spatial cut only as initial value (check if cluster is ellipse defined by cuts is done below): - double closestClusterDistance = sqrt(spatial_cuts_[detector].x() * spatial_cuts_[detector].x() + - spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + auto averageTimestamp = calculate_average_timestamp(seedTrack.get()); + seedTrack->setTimestamp(averageTimestamp); - double timeCut = std::max(time_cut_ref_track, time_cuts_[detector]); - LOG(DEBUG) << "Using timing cut of " << Units::display(timeCut, {"ns", "us", "s"}); + track->setTimestamp(averageTimestamp); + if(use_volume_scatterer_) { + track->setVolumeScatter(volume_radiation_length_); + } + seedTrack->setParticleMomentum(momentum_); + seedTrack->setMagneticField(magnetic_field_); + track->setParticleMomentum(momentum_); + track->setMagneticField(magnetic_field_); + + // Loop over each subsequent plane and look for a cluster within the timing cuts + size_t detector_nr = seed_detectors.size(); + // Get all detectors here to also include passive layers which might contribute to scattering + for(auto& detector : get_detectors()) { + if(detector->isAuxiliary()) { + continue; + } + auto detectorID = detector->getName(); + LOG(TRACE) << "added material budget for " << detectorID << " at z = " << detector->displacement().z(); + track->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); - auto neighbors = trees[detector].getAllElementsInTimeWindow(refTrack.timestamp(), timeCut); + if(std::find(seed_detectors.begin(), seed_detectors.end(), detector) != seed_detectors.end()) { + continue; + } - LOG(DEBUG) << "- found " << neighbors.size() << " neighbors within the correct time window"; + // register plane for seed track here to avoid doubly registering seed detectors + seedTrack->registerPlane(detectorID, detector->origin(), detector->materialBudget(), detector->toLocal()); - // Now look for the spatially closest cluster on the next plane - refTrack.fit(); + if(exclude_DUT_ && detector->isDUT()) { + LOG(DEBUG) << "Skipping DUT plane."; + continue; + } - PositionVector3D> interceptPoint = detector->getLocalIntercept(&refTrack); - double interceptX = interceptPoint.X(); - double interceptY = interceptPoint.Y(); + if(detector->isPassive()) { + LOG(DEBUG) << "Skipping passive plane."; + continue; + } - for(size_t ne = 0; ne < neighbors.size(); ne++) { - auto newCluster = neighbors[ne].get(); + // Determine whether a track can still be assembled given the number of current hits and the number of + // detectors to come. Reduces computing time. If this check fails, one can directly continue with the next seed + detector_nr++; + if(seedTrack->getNClusters() + (trees.size() - detector_nr + 1) < min_hits_on_track_) { + LOG(DEBUG) << "No chance to find a track - too few detectors left: " << seedTrack->getNClusters() << " + " + << trees.size() << " - " << detector_nr << " < " << min_hits_on_track_; + break; + } - // Calculate the distance to the previous plane's cluster/intercept - double distanceX = interceptX - newCluster->local().x(); - double distanceY = interceptY - newCluster->local().y(); - double distance = sqrt(distanceX * distanceX + distanceY * distanceY); + if(trees.count(detector) == 0) { + LOG(TRACE) << "Skipping detector " << detector->getName() << " as it has 0 clusters."; + continue; + } - // Check if newCluster lies within ellipse defined by spatial cuts around intercept, - // following this example: - // https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/ - // - // ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse, - // > 1: outside, - // < 1: inside - // Continue if outside of ellipse: + // Get all neighbors within the timing cut + LOG(DEBUG) << "Searching for neighboring cluster on device " << detector->getName(); + LOG(DEBUG) << "- reference time is " << Units::display(seedTrack->timestamp(), {"ns", "us", "s"}); + Cluster* closestCluster = nullptr; - double norm = (distanceX * distanceX) / (spatial_cuts_[detector].x() * spatial_cuts_[detector].x()) + - (distanceY * distanceY) / (spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + // Use spatial cut only as initial value (check if cluster is ellipse defined by cuts is done below): + double closestClusterDistance = sqrt(spatial_cuts_[detector].x() * spatial_cuts_[detector].x() + + spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); - if(norm > 1) { - LOG(DEBUG) << "Cluster outside the cuts. Normalized distance: " << norm; - continue; - } + double timeCut = std::max(time_cut_ref_track, time_cuts_[detector]); + LOG(DEBUG) << "Using timing cut of " << Units::display(timeCut, {"ns", "us", "s"}); - // If this is the closest keep it for now - if(distance < closestClusterDistance) { - closestClusterDistance = distance; - closestCluster = newCluster; - } - } + auto neighbors = trees[detector].getAllElementsInTimeWindow(seedTrack->timestamp(), timeCut); - if(closestCluster == nullptr) { - LOG(DEBUG) << "No cluster within spatial cut"; - continue; - } + LOG(DEBUG) << "- found " << neighbors.size() << " neighbors within the correct time window"; - // Add the cluster to the track - refTrack.addCluster(closestCluster); - track->addCluster(closestCluster); - averageTimestamp = calculate_average_timestamp(&refTrack); - refTrack.setTimestamp(averageTimestamp); - track->setTimestamp(averageTimestamp); - // At least 3 clusters, helix fit is possible - refTrack.setMagneticField(magnetic_field_); + // Now look for the spatially closest cluster on the next plane + seedTrack->fit(); - LOG(DEBUG) << "- added cluster to track"; - } + PositionVector3D> interceptPoint = detector->getLocalIntercept(seedTrack.get()); + double interceptX = interceptPoint.X(); + double interceptY = interceptPoint.Y(); - // check if track has required detector(s): - auto foundRequiredDetector = [this](Track* t) { - for(auto& requireDet : require_detectors_) { - if(!requireDet.empty() && !t->hasDetector(requireDet)) { - LOG(DEBUG) << "No cluster from required detector " << requireDet << " on the track."; - return false; - } + for(size_t ne = 0; ne < neighbors.size(); ne++) { + auto newCluster = neighbors[ne].get(); + + // Calculate the distance to the previous plane's cluster/intercept + double distanceX = interceptX - newCluster->local().x(); + double distanceY = interceptY - newCluster->local().y(); + double distance = sqrt(distanceX * distanceX + distanceY * distanceY); + + // Check if newCluster lies within ellipse defined by spatial cuts around intercept, + // following this example: + // https://www.geeksforgeeks.org/check-if-a-point-is-inside-outside-or-on-the-ellipse/ + // + // ellipse defined by: x^2/a^2 + y^2/b^2 = 1: on ellipse, + // > 1: outside, + // < 1: inside + // Continue if outside of ellipse: + + double norm = (distanceX * distanceX) / (spatial_cuts_[detector].x() * spatial_cuts_[detector].x()) + + (distanceY * distanceY) / (spatial_cuts_[detector].y() * spatial_cuts_[detector].y()); + + if(norm > 1) { + LOG(DEBUG) << "Cluster outside the cuts. Normalized distance: " << norm; + continue; + } + + // If this is the closest keep it for now + if(distance < closestClusterDistance) { + closestClusterDistance = distance; + closestCluster = newCluster; } - return true; - }; - if(!foundRequiredDetector(track.get())) { - continue; } - // Now should have a track with one cluster from each plane - if(track->getNClusters() < min_hits_on_track_) { - LOG(DEBUG) << "Not enough clusters on the track, found " << track->getNClusters() << " but " - << min_hits_on_track_ << " required."; + if(closestCluster == nullptr) { + LOG(DEBUG) << "No cluster within spatial cut"; continue; } - // Fit the track - track->fit(); + // Add the cluster to the track + seedTrack->addCluster(closestCluster); + track->addCluster(closestCluster); + averageTimestamp = calculate_average_timestamp(seedTrack.get()); + seedTrack->setTimestamp(averageTimestamp); + track->setTimestamp(averageTimestamp); - if(reject_by_ROI_ && track->isFitted()) { - // check if the track is within ROI for all detectors - auto ds = get_regular_detectors(!exclude_DUT_); - auto out_of_roi = - std::find_if(ds.begin(), ds.end(), [track](const auto& d) { return !d->isWithinROI(track.get()); }); - if(out_of_roi != ds.end()) { - LOG(DEBUG) << "Rejecting track outside of ROI of detector " << out_of_roi->get()->getName(); - continue; + LOG(DEBUG) << "- added cluster to track"; + } + + // check if track has required detector(s): + auto foundRequiredDetector = [this](Track* t) { + for(auto& requireDet : require_detectors_) { + if(!requireDet.empty() && !t->hasDetector(requireDet)) { + LOG(DEBUG) << "No cluster from required detector " << requireDet << " on the track."; + return false; } } - // save the track - if(track->isFitted()) { - tracks.push_back(track); - } else { - LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting"; + return true; + }; + if(!foundRequiredDetector(track.get())) { + continue; + } + + // Now should have a track with one cluster from each plane + if(track->getNClusters() < min_hits_on_track_) { + LOG(DEBUG) << "Not enough clusters on the track, found " << track->getNClusters() << " but " + << min_hits_on_track_ << " required."; + continue; + } + + // Fit the track + track->fit(); + + if(reject_by_ROI_ && track->isFitted()) { + // check if the track is within ROI for all detectors + auto ds = get_regular_detectors(!exclude_DUT_); + auto out_of_roi = + std::find_if(ds.begin(), ds.end(), [track](const auto& d) { return !d->isWithinROI(track.get()); }); + if(out_of_roi != ds.end()) { + LOG(DEBUG) << "Rejecting track outside of ROI of detector " << out_of_roi->get()->getName(); continue; } + } + // save the track + if(track->isFitted()) { + tracks.push_back(track); + } else { + LOG_N(WARNING, 100) << "Rejected a track due to failure in fitting"; + continue; + } - if(timestamp_from_.empty()) { - // Improve the track timestamp by taking the average of all planes - auto timestamp = calculate_average_timestamp(track.get()); - track->setTimestamp(timestamp); - LOG(DEBUG) << "Using average cluster timestamp of " << Units::display(timestamp, "us") - << " as track timestamp."; - } else { - // use timestamp of required detector: - double track_timestamp = track->getClusterFromDetector(timestamp_from_)->timestamp(); - LOG(DEBUG) << "Using timestamp of detector " << timestamp_from_ - << " as track timestamp: " << Units::display(track_timestamp, "us"); - track->setTimestamp(track_timestamp); - } + if(timestamp_from_.empty()) { + // Improve the track timestamp by taking the average of all planes + auto timestamp = calculate_average_timestamp(track.get()); + track->setTimestamp(timestamp); + LOG(DEBUG) << "Using average cluster timestamp of " << Units::display(timestamp, "us") << " as track timestamp."; + } else { + // use timestamp of required detector: + double track_timestamp = track->getClusterFromDetector(timestamp_from_)->timestamp(); + LOG(DEBUG) << "Using timestamp of detector " << timestamp_from_ + << " as track timestamp: " << Units::display(track_timestamp, "us"); + track->setTimestamp(track_timestamp); } } -- GitLab From 5cd4611b4ebe94b1debf40b50e18b0b83095780c Mon Sep 17 00:00:00 2001 From: David Bacher Date: Thu, 14 Jul 2022 13:37:55 +0200 Subject: [PATCH 12/14] HelixTrack: Remove straight lines from helix model --- src/objects/HelixTrack.cpp | 184 ++++++++++--------------------------- src/objects/HelixTrack.hpp | 2 - 2 files changed, 50 insertions(+), 136 deletions(-) diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp index 58d8964b..62b06768 100644 --- a/src/objects/HelixTrack.cpp +++ b/src/objects/HelixTrack.cpp @@ -45,9 +45,8 @@ ROOT::Math::XYZVector HelixTrack::getDirection(const double& z) const { void HelixTrack::calculate_transformations() { // compute transformation to helix coordinate system with the direction of the magnetic field as the z-direction - // If there is no magnetic field, the transformation will simply result in identity XYZVector globalZ = XYZVector(0, 0, 1); - XYZVector helixZ = magneticField_.R() != 0 ? magneticField_.unit() : globalZ; + XYZVector helixZ = magneticField_.unit(); XYZVector helixY = helixZ.Theta() != 0 ? helixZ.Cross(globalZ).unit() : XYZVector(0, 1, 0); XYZVector helixX = helixY.Cross(helixZ).unit(); @@ -57,14 +56,9 @@ void HelixTrack::calculate_transformations() { void HelixTrack::calculate_chi2() { - // Straightline: Each hit provides two measurements and there are four fitting parameters. - if(magneticField_.R() == 0) { - ndof_ = (track_clusters_.size() - 2) * 2; - } - // Helix: Each hit provides two measurements and there are five fitting paramters. - else { - ndof_ = 2 * track_clusters_.size() - 5; - } + // Each hit provides two measurements and there are five fitting paramters. + ndof_ = 2 * track_clusters_.size() - 5; + chi2_ = 0.; chi2ndof_ = 0.; @@ -117,8 +111,16 @@ double HelixTrack::operator()(const double* parameters) { return chi2_; } -void HelixTrack::fit_helix() { - LOG(TRACE) << "Starting helix fit."; +void HelixTrack::fit() { + + if(magneticField_ == XYZVector()) { + throw TrackError(typeid(HelixTrack), " can only be used with magnetic fields"); + } + + isFitted_ = false; + calculate_transformations(); + + LOG(DEBUG) << "Starting helix fit."; // Fit points to helix with fast helix fit described in: // https://www.physi.uni-heidelberg.de/Publications/DiplomaKiehn.pdf @@ -139,7 +141,6 @@ void HelixTrack::fit_helix() { } auto position = toHelix_ * cluster->global(); - LOG(TRACE) << position; long double x = position.x(); long double y = position.y(); long double r2 = x * x + y * y; @@ -179,7 +180,7 @@ void HelixTrack::fit_helix() { curvature_ = static_cast(2. * kappa / sqrt(1. - 4. * delta * kappa)); dca_ = static_cast(2. * delta / (1. + sqrt(1. - 4. * delta * kappa))); - LOG(TRACE) << "Circle fit finished. Radius: " << Units::display(1 / curvature_, {"m", "km"}) + LOG(DEBUG) << "Circle fit finished. Radius: " << Units::display(1 / curvature_, {"m", "km"}) << ", DCA: " << Units::display(dca_, {"mm", "um"}) << ", Phi: " << Units::display(phi_, "deg"); // 2. Calculate arc length in transverse plane @@ -213,78 +214,7 @@ void HelixTrack::fit_helix() { z0_ = res[0]; tanLambda_ = res[1]; - LOG(TRACE) << "Longitudinal fit finished. Z0: " << Units::display(z0_, {"mm", "cm"}) << ", tan(lambda): " << tanLambda_; -} - -void HelixTrack::fit_straightLine() { - LOG(TRACE) << "No Magnetic Field. Starting straightline fit."; - Matrix4d mat(Matrix4d::Zero()); - Vector4d vec(Vector4d::Zero()); - - for(auto& cl : track_clusters_) { - auto* cluster = cl.get(); - if(cluster == nullptr) { - throw MissingReferenceException(typeid(*this), typeid(Cluster)); - } - - double x = cluster->global().x(); - double y = cluster->global().y(); - double z = cluster->global().z(); - double ex2 = cluster->errorX() * cluster->errorX(); - double ey2 = cluster->errorY() * cluster->errorY(); - - Vector2d meas(x, y); - Matrix2d V; - V << ex2, 0., 0., ey2; - Matrix C; - C << 1., z, 0., 0., 0., 0., 1., z; - - vec += C.transpose() * V.inverse() * meas; - mat += C.transpose() * V.inverse() * C; - } - - // Check for singularities. - if(fabs(mat.determinant()) < std::numeric_limits::epsilon()) { - throw TrackFitError(typeid(this), "Martix inversion in straight line fit failed"); - } - // Get the StraightLineTrack parameters - Eigen::Vector4d res = mat.inverse() * vec; - - double bx = res[0]; - double by = res[2]; - double mx = res[1]; - double my = res[3]; - - LOG(TRACE) << "bx = " << bx << ", by = " << by << ", mx = " << mx << ", my = " << my; - - // cover case without any slope in x or y direction - if(mx < std::numeric_limits::epsilon() && my < std::numeric_limits::epsilon()) { - tanLambda_ = std::numeric_limits::max(); - phi_ = atan2(bx, -by); - dca_ = sin(phi_) * bx - cos(phi_) * by; - z0_ = 0.; - } else { - tanLambda_ = 1. / sqrt(mx * mx + my * my); - phi_ = atan2(my, mx); - dca_ = tanLambda_ * (bx * my - by * mx); - z0_ = -(bx * mx + by * my) * tanLambda_ * tanLambda_; - } - curvature_ = 0.; - LOG(TRACE) << "Straightline fit finished. DCA: " << Units::display(dca_, {"m", "cm", "mm", "um"}) - << ", Phi: " << Units::display(phi_, "deg") << ", z0: " << Units::display(z0_, {"m", "cm", "mm", "um"}) - << ", tan(lambda): " << tanLambda_; -} - -void HelixTrack::fit() { - - isFitted_ = false; - calculate_transformations(); - - if(magneticField_.R() == 0.) { - fit_straightLine(); - } else { - fit_helix(); - } + LOG(DEBUG) << "Longitudinal fit finished. Z0: " << Units::display(z0_, {"mm", "cm"}) << ", tan(lambda): " << tanLambda_; // Calculate the chi2 this->calculate_chi2(); @@ -299,11 +229,6 @@ ROOT::Math::XYZPoint HelixTrack::getIntercept(double z) const { } double HelixTrack::getArcLength(double x, double y) const { - // Straightline - if(curvature_ == 0.) { - return cos(phi_) * x + sin(phi_) * y; - } - // Helix // Transform to cooridnate system with origin in center of track circle Vector2d posOrigin; posOrigin << x, y; @@ -326,32 +251,25 @@ ROOT::Math::XYZPoint HelixTrack::getIntercept(XYZPoint origin, XYZVector normal) double sinLambda = sqrt(1. - cosLambda * cosLambda); XYZVector distance, trackDirection; XYZPoint position; - if(curvature_ == 0.) { - trackDirection = XYZVector(cosLambda * cos(phi_), cosLambda * sin(phi_), sinLambda); - // Point of closest approach - XYZPoint pca(dca_ * sin(phi_), -dca_ * cos(phi_), z0_); - distance = pca - originHelix; - double pathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); - position = pca + pathLength * trackDirection; - } else { - double a = getArcLength(originHelix.x(), originHelix.y()); - for(unsigned int nIter = 0; nIter < 10; nIter++) { - const double dPhi = a * curvature_; - const double cosPhi = cos(dPhi - phi_); - const double sinPhi = sin(dPhi - phi_); - trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); - position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, - -cos(phi_) * (dca_ + 1. / curvature_) + cosPhi * 1. / curvature_, - z0_ + tanLambda_ * a); - distance = position - originHelix; - const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); - if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um - a += corrPathLength * cosLambda; - } else { - break; - } + // initial estimate of arc length + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(dPhi - phi_); + const double sinPhi = sin(dPhi - phi_); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) + cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; } } + return toGlobal_ * position; } @@ -364,27 +282,25 @@ ROOT::Math::XYZVector HelixTrack::getDirection(XYZPoint origin, XYZVector normal double sinLambda = sqrt(1. - cosLambda * cosLambda); XYZVector distance, trackDirection; XYZPoint position; - if(curvature_ == 0.) { - trackDirection = XYZVector(cosLambda * cos(phi_), cosLambda * sin(phi_), sinLambda); - } else { - double a = getArcLength(originHelix.x(), originHelix.y()); - for(unsigned int nIter = 0; nIter < 10; nIter++) { - const double dPhi = a * curvature_; - const double cosPhi = cos(phi_ + dPhi); - const double sinPhi = sin(phi_ + dPhi); - trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); - position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, - -cos(phi_) * (dca_ + 1. / curvature_) - cosPhi * 1. / curvature_, - z0_ + tanLambda_ * a); - distance = position - originHelix; - const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); - if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um - a += corrPathLength * cosLambda; - } else { - break; - } + // initial estimate of arc length + double a = getArcLength(originHelix.x(), originHelix.y()); + for(unsigned int nIter = 0; nIter < 10; nIter++) { + const double dPhi = a * curvature_; + const double cosPhi = cos(phi_ + dPhi); + const double sinPhi = sin(phi_ + dPhi); + trackDirection = XYZVector(cosLambda * cosPhi, cosLambda * sinPhi, sinLambda); + position = XYZPoint(sin(phi_) * (dca_ + 1. / curvature_) + sinPhi * 1. / curvature_, + -cos(phi_) * (dca_ + 1. / curvature_) - cosPhi * 1. / curvature_, + z0_ + tanLambda_ * a); + distance = position - originHelix; + const double corrPathLength = -distance.Dot(normalHelix) / trackDirection.Dot(normalHelix); + if(fabs(corrPathLength) > 1e-4) { // Correction is more than 0.1 um + a += corrPathLength * cosLambda; + } else { + break; } } + return toGlobal_ * trackDirection; } diff --git a/src/objects/HelixTrack.hpp b/src/objects/HelixTrack.hpp index f7f84c5f..60581b9b 100644 --- a/src/objects/HelixTrack.hpp +++ b/src/objects/HelixTrack.hpp @@ -106,8 +106,6 @@ namespace corryvreckan { void calculate_transformations(); - void fit_straightLine(); - void fit_helix(); // Member variables (Karimaeki parametrization) double curvature_; double phi_; -- GitLab From ff71b81a80920a6f3dc34c0b2b74200e1e071bfe Mon Sep 17 00:00:00 2001 From: David Bacher Date: Thu, 14 Jul 2022 13:44:29 +0200 Subject: [PATCH 13/14] Tracking4D: throw error if helix is used without magnetic field --- src/modules/Tracking4D/Tracking4D.cpp | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index 0f4e82db..74a9ce25 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -71,10 +71,15 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("unique_cluster_usage"); magnetic_field_ = config_.get("magnetic_field"); - // straightline tracks with magnetic field are unphysical + // the folowing combinations of track model and magnetic field are not supported: + // straightline with magnetic field + // helix without magnetic field if(track_model_ == "straightline" && magnetic_field_ != ROOT::Math::XYZVector()) { throw InvalidCombinationError( - config_, {"track_model", "magnetic_field"}, "Straightlines do not support magnetic fields."); + config_, {"track_model", "magnetic_field"}, track_model_ + " does not support magnetic fields."); + } else if(track_model_ == "helix" && magnetic_field_ == ROOT::Math::XYZVector()) { + throw InvalidCombinationError( + config_, {"track_model", "magnetic_field"}, track_model_ + " can only be used with magnetic fields."); } // print a warning if volumeScatterer are used as this causes fit failures -- GitLab From 0d4e9d0d444840094806023c4cbfd4c603917b40 Mon Sep 17 00:00:00 2001 From: David Bacher Date: Tue, 19 Jul 2022 09:45:47 +0200 Subject: [PATCH 14/14] Tracking4D: add check on radius during track reconstruction --- src/modules/Tracking4D/Tracking4D.cpp | 44 +++++++++++++++++++++++++++ src/modules/Tracking4D/Tracking4D.h | 5 +++ src/objects/HelixTrack.cpp | 7 +++++ src/objects/HelixTrack.hpp | 7 +++++ src/objects/Track.cpp | 4 +++ src/objects/Track.hpp | 7 +++++ 6 files changed, 74 insertions(+) diff --git a/src/modules/Tracking4D/Tracking4D.cpp b/src/modules/Tracking4D/Tracking4D.cpp index 74a9ce25..bd99a2cf 100644 --- a/src/modules/Tracking4D/Tracking4D.cpp +++ b/src/modules/Tracking4D/Tracking4D.cpp @@ -36,6 +36,8 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("reject_by_roi", false); config_.setDefault("unique_cluster_usage", false); config_.setDefault("magnetic_field", ROOT::Math::XYZVector()); + config_.setDefault("radius_tol_min", 1.0); + config_.setDefault("radius_tol_rel", 0.1); if(config_.count({"time_cut_rel", "time_cut_abs"}) == 0) { config_.setDefault("time_cut_rel", 3.0); @@ -70,6 +72,16 @@ Tracking4D::Tracking4D(Configuration& config, std::vector("reject_by_roi"); unique_cluster_usage_ = config_.get("unique_cluster_usage"); magnetic_field_ = config_.get("magnetic_field"); + radius_tol_min_ = config_.get("radius_tol_min"); + radius_tol_rel_ = config_.get("radius_tol_rel"); + + expected_radius_ = 0.; + if(magnetic_field_ != ROOT::Math::XYZVector()) { + expected_radius_ = momentum_ * + sqrt(magnetic_field_.x() * magnetic_field_.x() + magnetic_field_.y() * magnetic_field_.y()) / + (magnetic_field_.R() * magnetic_field_.R() * 0.3); // R[mm] = pT[MeV]/(B[T]*c)*1000 + } + radius_tol_ = std::max(radius_tol_min_, expected_radius_ * radius_tol_rel_); // the folowing combinations of track model and magnetic field are not supported: // straightline with magnetic field @@ -110,6 +122,12 @@ void Tracking4D::initialize() { trackAngleX = new TH1F("trackAngleX", title.c_str(), 2000, -0.01, 0.01); title = "Track angle Y;angle_{y} [rad];events"; trackAngleY = new TH1F("trackAngleY", title.c_str(), 2000, -0.01, 0.01); + title = "Track radius;radius [m];events"; + trackRadius = new TH1F("trackRadius", + title.c_str(), + 500, + std::max(0., (expected_radius_ - radius_tol_) / 1000), + (expected_radius_ + radius_tol_) / 1000); title = "Track time within event;track time - event start;events"; trackTime = new TH1F("trackTime", title.c_str(), 1000, 0, 460.8); title = "Track time with respect to first trigger;track time - trigger;events"; @@ -474,6 +492,19 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { // Now look for the spatially closest cluster on the next plane seedTrack->fit(); + // check whether the seed track has a radius within the limits of the expectation from momentum and magnetic + // field. Reduces computing time + if(magnetic_field_ != ROOT::Math::XYZVector()) { + auto fittedRadius = seedTrack->getRadius(); + if(std::fabs(fittedRadius - expected_radius_) > radius_tol_) { + LOG(DEBUG) << "Radius of seed track outside of tolerance. Looking at next seed. Expected radius: " + << Units::display(expected_radius_, {"mm", "m", "km", "um"}) + << ", tolerance: " << Units::display(radius_tol_, {"mm", "m", "km", "um"}) + << ", fitted radius: " << Units::display(fittedRadius, {"mm", "m", "km", "um"}); + break; + } + } + PositionVector3D> interceptPoint = detector->getLocalIntercept(seedTrack.get()); double interceptX = interceptPoint.X(); double interceptY = interceptPoint.Y(); @@ -559,6 +590,18 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { continue; } } + // check whether the track has a radius within the limits of the expectation from momentum and magnetic field. + // Reduces computing time + if(magnetic_field_ != ROOT::Math::XYZVector() && track->isFitted() && track_model_ == "helix") { + auto fittedRadius = track->getRadius(); + if(std::fabs(fittedRadius - expected_radius_) > radius_tol_) { + LOG(DEBUG) << "Rejecting track with radius outside of tolerance. Expected radius: " + << Units::display(expected_radius_, {"mm", "m", "km", "um"}) + << ", tolerance: " << Units::display(radius_tol_, {"mm", "m", "km", "um"}) + << ", fitted radius: " << Units::display(fittedRadius, {"mm", "m", "km", "um"}); + continue; + } + } // save the track if(track->isFitted()) { tracks.push_back(track); @@ -639,6 +682,7 @@ StatusCode Tracking4D::run(const std::shared_ptr& clipboard) { trackAngleX->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).X())); trackAngleY->Fill(atan(track->getDirection(track->getClusters().front()->detectorID()).Y())); } + trackRadius->Fill(track->getRadius() / 1000); // Make residuals auto trackClusters = track->getClusters(); for(auto& trackCluster : trackClusters) { diff --git a/src/modules/Tracking4D/Tracking4D.h b/src/modules/Tracking4D/Tracking4D.h index d13a8c5d..b2417ed6 100644 --- a/src/modules/Tracking4D/Tracking4D.h +++ b/src/modules/Tracking4D/Tracking4D.h @@ -46,6 +46,7 @@ namespace corryvreckan { TH1F* tracksPerEvent; TH1F* trackAngleX; TH1F* trackAngleY; + TH1F* trackRadius; TH1F* tracksVsTime; std::map residualsX_local; std::map residualsXwidth1_local; @@ -95,6 +96,10 @@ namespace corryvreckan { std::string timestamp_from_; std::string track_model_; ROOT::Math::XYZVector magnetic_field_; + double expected_radius_; + double radius_tol_; + double radius_tol_min_; + double radius_tol_rel_; // Function to calculate the weighted average timestamp from the clusters of a track double calculate_average_timestamp(const Track* track); diff --git a/src/objects/HelixTrack.cpp b/src/objects/HelixTrack.cpp index 62b06768..69ebc8b5 100644 --- a/src/objects/HelixTrack.cpp +++ b/src/objects/HelixTrack.cpp @@ -23,6 +23,13 @@ ROOT::Math::XYPoint HelixTrack::getKinkAt(const std::string&) const { return ROOT::Math::XYPoint(0, 0); } +double HelixTrack::getRadius() const { + if(!isFitted_) { + throw TrackError(typeid(HelixTrack), " has no defined radius before fitting"); + } + return fabs(1 / curvature_); +} + ROOT::Math::XYZPoint HelixTrack::getState(const std::string& detectorID) const { const Plane* plane = get_plane(detectorID); XYZVector planeU, planeV, planeN; diff --git a/src/objects/HelixTrack.hpp b/src/objects/HelixTrack.hpp index 60581b9b..06c40bba 100644 --- a/src/objects/HelixTrack.hpp +++ b/src/objects/HelixTrack.hpp @@ -86,6 +86,13 @@ namespace corryvreckan { */ ROOT::Math::XYPoint getKinkAt(const std::string& detectorID) const override; + /** + * @brief Get the radius in the transverse plane at a given detector layer. + * @param detectorID Detector ID at which the kink should be evaluated + * @return radius in transverse plane at a given detector + */ + double getRadius() const override; + /** * @brief Get the arc length in the transverse plane of the helix * @param x x-coordinate in transverse plane diff --git a/src/objects/Track.cpp b/src/objects/Track.cpp index 31b352b3..2e618c04 100644 --- a/src/objects/Track.cpp +++ b/src/objects/Track.cpp @@ -231,6 +231,10 @@ double Track::getMaterialBudget(const std::string& detectorID) const { return budget; } +double Track::getRadius() const { + return 0.; +} + void Track::registerPlane(const std::string& name, ROOT::Math::XYZPoint origin, double x0, Transform3D g2l) { Plane p(name, origin, x0, g2l); auto pl = diff --git a/src/objects/Track.hpp b/src/objects/Track.hpp index 2fcfcbde..8e738c2a 100644 --- a/src/objects/Track.hpp +++ b/src/objects/Track.hpp @@ -229,6 +229,13 @@ namespace corryvreckan { */ virtual ROOT::Math::XYPoint getKinkAt(const std::string& detectorID) const = 0; + /** + * @brief Get the radius in the transverse plane at a given detector layer. + * @param detectorID Detector ID at which the kink should be evaluated + * @return radius in transverse plane at a given detector + */ + virtual double getRadius() const; + /** * @brief Get the materialBudget of a detector layer * @param detectorID -- GitLab